Não foi possível enviar o arquivo. Será algum problema com as permissões?
Essa é uma revisão anterior do documento!
Ridículas - dicas curtas sobre R
Ridículas é a página do LEG dedicada à fornecer dicas curtas sobre R, e.g. condução de análises, operação com dados e confecção de gráficos. As dicas estão organizadas pelo título, seguido de descrição, palavras-chave e CMR (código mínimo reproduzível). Se você deseja contribuir com a nossa página de Ridículas, envie e-mail para walmes@ufpr.br.
Regressão na análise de variância
No exemplo abaixo fazemos a análise de um experimento fatorial. São combinados três níveis de um fator qualitativo (cultivar de sorgo) com 6 níveis de um fator quantitativo (dose de nitrogênio). A resposta avaliada nas parcelas é o índice agronômico. O experimento foi instalado em blocos. A interação entre as fontes de variação foi significativa. O desdobramento dessa interação foi feita nos dois sentidos, ou seja, estudo do fator dose fixando nível do fator cultivar e estudo do fator cultivar fixando o nível do fator dose. Foi feito a ajuste de modelos de regressão polinomial para o estudo do fator dose. São ilustrados passos para teste de hipótese sobre os efeitos linear, quadrático, etc, e falta de ajuste do modelo ajustado, além de como obter as estimativas dos termos do modelo e o coeficiente de determinação (R²).
palavras-chave: #regressão, #polinômio, #anova.
#------------------------------------------------------------------------------------------
# por Walmes
#------------------------------------------------------------------------------------------
# dados
sorgo <- read.table("http://www.leg.ufpr.br/~walmes/docs/anovareg.txt", header=TRUE)
sorgo <- transform(sorgo, bloco=factor(bloco), cultivar=factor(cultivar))
str(sorgo)
#
#------------------------------------------------------------------------------------------
# gráficos exploratórios
require(lattice)
xyplot(indice~dose|cultivar, groups=bloco, data=sorgo,
jitter.x=TRUE, type=c("p","l"), layout=c(3,1))
xyplot(indice~dose, groups=cultivar, data=sorgo, jitter.x=TRUE, type=c("p","a"))
#
#------------------------------------------------------------------------------------------
# análise de variância do modelo de fatores
m0 <- aov(indice~bloco+cultivar*ordered(dose), data=sorgo)
summary(m0)
#
#------------------------------------------------------------------------------------------
# checagem
par(mfrow=c(2,2))
plot(m0)
layout(1)
#
#------------------------------------------------------------------------------------------
# desdobrando as somas de quadrados de doses dentro de cultivar
# dicas: forneça para ’by’ o número de níveis de cultivar (3)
# forneça para ’length.out’ os graus de liberdade de dose (6-1)
m1 <- aov(indice~bloco+cultivar/ordered(dose), data=sorgo)
summary(m1)
coef(m1)
summary(m1, split=list("cultivar:ordered(dose)"=list(
"Ag-1002"=seq(1, by=3, length.out=5),
"BR-300"=seq(2, by=3, length.out=5),
"Pioneer-B815"=seq(3, by=3, length.out=5)
)))
#
#------------------------------------------------------------------------------------------
# desdobrando somas de quadrados de cultivar dentro das doses
# dicas: forneça para ’by’ o número de níveis de dose (6)
# forneça para ’length.out’ os graus de liberdade de cultivar (3-1)
m2 <- aov(indice~bloco+ordered(dose)/cultivar, data=sorgo)
coef(m2)
summary(m2, split=list("ordered(dose):cultivar"=list(
"N.0"=seq(1, by=6, length.out=2),
"N.60"=seq(2, by=6, length.out=2),
"N.120"=seq(3, by=6, length.out=2),
"N.180"=seq(4, by=6, length.out=2),
"N.240"=seq(5, by=6, length.out=2),
"N.300"=seq(6, by=6, length.out=2)
)))
#
#------------------------------------------------------------------------------------------
# desdobrando efeitos dos graus polinômio dentro de dose dentro de cultivar
# lof é falta de ajuste (lack of fit)
summary(m1, split=list("cultivar:ordered(dose)"=list(
"Ag-1002.L"=1,
"Ag-1002.Q"=4,
"Ag-1002.C"=7,
"Ag-1002.lof"=c(10,13),
"BR-300.L"=2,
"BR-300.Q"=5,
"BR-300.C"=8,
"BR-300.lof"=c(11,14),
"Pioneer-B815.L"=3,
"Pioneer-B815.Q"=6,
"Pioneer-B815.C"=9,
"Pioneer-B815.lof"=c(12,15)
)))
#
#------------------------------------------------------------------------------------------
# obter as equações de regressão e R^2 para os modelos linear, quadrático e cúbico
# dica: usar contraste tipo soma zero para blocos para se anularem na fórmula
# e remover o intercepto especificando o ’-1’, trocar a ordem dos termos no modelo
# linear (estimativas corretas mas erros padrões e p-valores precisam de correção)
m3 <- aov(indice~-1+cultivar/dose+bloco, data=sorgo,
contrast=list(bloco=contr.sum))
summary.lm(m3)
#
#------------------------------------------------------------------------------------------
# quadrático (estimativas corretas mas erros padrões e p-valores precisam de correção)
m4 <- aov(indice~-1+cultivar/(dose+I(dose^2))+bloco, data=sorgo,
contrast=list(bloco=contr.sum))
summary.lm(m4)
#
#------------------------------------------------------------------------------------------
# cúbico (estimativas corretas mas erros padrões e p-valores precisam de correção)
m5 <- aov(indice~-1+cultivar/(dose+I(dose^2)+I(dose^3))+bloco, data=sorgo,
contrast=list(bloco=contr.sum))
summary.lm(m5)
#
#------------------------------------------------------------------------------------------
# calcular os R^2
sapply(c(linear=1, quadrático=2, cúbico=3),
function(degree){
sapply(levels(sorgo$cultivar),
function(i){
da <- with(subset(sorgo, cultivar==i),
aggregate(indice, list(dose=dose), mean))
summary(lm(x~poly(dose, degree, raw=TRUE), da))$r.squared
})})
#
#------------------------------------------------------------------------------------------
Experimento com dois fatores de efeito aditivo e perda de muitas parcelas
Apresenta a análise de um experimento de dois fatores qualitativos com a perda de muitas parcelas. Nesse caso, a perda de parcelas não compromete a estimabilidade dos efeitos dos níveis dos fatores, apenas implica que eles sejam estimados com precisões diferentes. Além disso, a ortogonalidade entre os efeitos é perdida e a interpretação dos testes de hipótese pela análise de variância requer cuidados. Finalmente, é ilustrado a obtenção das estimativas das médias marginais para os níveis dos fatores (conhecidas por lsmeans na documentação de outro aplicativo) e realizadas as comparações múltiplas dessas médias com correção para o nível nominal de significância.
palavras-chave: #parcela_perdida, #desbalanceamento, #médias_ajustadas, #aditivo.
#------------------------------------------------------------------------------------------
# por Walmes
#------------------------------------------------------------------------------------------
# dados
da <- expand.grid(rept=1:5, ep=factor(1:5), tr=factor(1:4))
da$y <- c(58.4, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
68.4, NA, NA, NA, NA, 258.8, 265.6, NA, NA, NA, NA, NA, 250, NA, 278.8,
268.8, NA, NA, NA, 309.6, NA, NA, NA, NA, NA, NA, NA, NA, NA, 254, 598.8,
NA, NA, NA, NA, 250, 399.6, 260, NA, NA, NA, 288.4, NA, NA, NA, 397.2, NA,
NA, 337.6, NA, 415.2, NA, 450.8, NA, NA, NA, NA, 393.2, NA, NA, NA, NA,
NA, NA, NA, 380.4, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 634, 417.2,
NA, NA, NA, NA, NA)
#------------------------------------------------------------------------------------------
# ajuste do modelo aditivo com teste F marginal
m0 <- lm(y~ep+tr, data=da)
drop1(m0, test="F")
#------------------------------------------------------------------------------------------
# análise gráfica dos resíduos
par(mfrow=c(2,2))
plot(m0)
layout(1)
#------------------------------------------------------------------------------------------
# estimativas dos efeitos (solução) sob a restrição do R
summary(m0)
#------------------------------------------------------------------------------------------
# obtenção das médias ajustadas dos níveis de tratamento
require(contrast)
lapply(levels(da$tr),
function(i){
contrast(m0, type="average", list(tr=i, ep=levels(da$ep)))
}
)
#------------------------------------------------------------------------------------------
# comparação múltipla de efeitos
require(multcomp)
summary(glht(m0, linfct=mcp(tr="Tukey")))
#------------------------------------------------------------------------------------------
Experimento em parcelas subdivididas com resposta do tipo binomial
Temporariamente sem descrição.
palavras-chave: #binomial, #subdividida, #verossimilhança, #experimento, #glm.
#------------------------------------------------------------------------------------------
# por Walmes
#------------------------------------------------------------------------------------------
# pacote que contem a função glmmPQL. pode-se usar a lme4::glmer()
require(MASS)
#------------------------------------------------------------------------------------------
# simulandos com a estrutura de um modelo de parcelas subdivididas e resposta binomial
# melhor forma de entender como são gerados os valores observados num experimento
nA <- 4; # número de níveis de um fator A de efeito fixo
nB <- 8; # número de níveis de um fator B de efeito aleatório
nC <- 2; # número de níveis de um fator C de efeito fixo
da <- expand.grid(A=factor(1:nA), B=factor(1:nB), C=factor(1:nC))
da$parcela <- with(da, factor(paste(A,B,sep="-"))) # níveis de parcela aleatório
XA <- model.matrix(~A-1, data=da) # matriz dos efeitos fixos para A
XB <- model.matrix(~B-1, data=da) # matriz dos efeitos aleatórios de B
XC <- model.matrix(~C-1, data=da) # matriz dos efeitos fixos para C
Xp <- model.matrix(~parcela-1, data=da) # matriz dos efeitos aleatórios de parcela
bA <- c(-1,0,1,1) # vetor de efeitos fixos A
bB <- rnorm(nlevels(da$B),0,1) # vetor de efeitos aleatórios B
bC <- c(0,-1) # vetor de efeitos fixos C
bp <- rnorm(nlevels(da$parcela),0,1) # vetor de efeitos aleatórios parcela
da$eta <- XA%*%bA+XB%*%bB+XC%*%bC+Xp%*%bp # vetor de valores preditos
summary(exp(da$eta)/(1+exp(da$eta))) # vetor das probabilidades preditas
da$yobs <- rbinom(nrow(da),
prob=exp(da$eta)/(1+exp(da$eta)),
size=100) # vetor de dados observados
summary(da$yobs)
str(da)
#------------------------------------------------------------------------------------------
# estimação de parâmetros via máxima verossimilhança (usando Penalized Quasi Likelihood)
m0 <- glmmPQL(cbind(yobs, 100-yobs)~A*C, random=~1|B/parcela, data=da, family=binomial)
summary(m0)
#------------------------------------------------------------------------------------------
Duas ou mais funções em um mesmo gráfico
Temporariamente sem descrição.
palavras-chave: #gráfico, #função, #curva.
#------------------------------------------------------------------------------------------ # por Walmes #------------------------------------------------------------------------------------------ # cria os vetores usando as funções abaixo e plota 2 num mesmo gráfico x <- seq(0, 10, by=0.1) y1 <- 4-0.3*x+0.02*x^2 y2 <- 3+exp(-0.9*x) matplot(x, cbind(y1, y2), lty=1, col=1:2, type="l") #------------------------------------------------------------------------------------------
Representação de ajuste de curvas de regressão no mesmo gráfico
Faremos o ajuste de três curvas de regressão à dados de índice agronômico de sorgo em função da dose de nitrogênio aplicada. Os dados estão classificados para três cultivares de sorgo. Apesar de se tratar de um experimento em blocos, isso será ignorado no exemplo que tem o objetivo de apresentar os comandos para ajustar curvas e representar o ajuste no mesmo gráfico. palavras-chave: #regressão, #lapply, #curvas, #ajuste, #lines, #subset.
#------------------------------------------------------------------------------------------
# por Walmes
#------------------------------------------------------------------------------------------
# dados.
sorgo <- read.table("http://www.leg.ufpr.br/~walmes/docs/anovareg.txt", header=TRUE)
sorgo <- transform(sorgo, bloco=factor(bloco), cultivar=factor(cultivar))
str(sorgo)
#------------------------------------------------------------------------------------------
# usaremos esse banco de dados para ajustarmos curvas de regressão separadas por cultivar.
# iremos ignorar a presença de blocos assumindo amostragem aleatória simples dos dados.
#------------------------------------------------------------------------------------------
# gráfico de dispersão dos dados.
levels(sorgo$cultivar)
par(mfrow=c(1,3))
with(sorgo, sapply(levels(sorgo$cultivar),
function(x){
plot(indice~dose, subset(sorgo, cultivar==x), main=x)
}
)
)
#------------------------------------------------------------------------------------------
# ajustaremos uma curva para cada cultivar. pode-se usar a lapply o tempo todo para isso.
levels(sorgo$cultivar)
mAG <- lm(indice~dose+I(dose^2), data=subset(sorgo, cultivar==levels(sorgo$cultivar)[1]))
mBR <- lm(indice~dose+I(dose^2), data=subset(sorgo, cultivar==levels(sorgo$cultivar)[2]))
mPi <- lm(indice~dose+I(dose^2), data=subset(sorgo, cultivar==levels(sorgo$cultivar)[3]))
lapply(list(mAG, mBR, mPi), summary)
lapply(list(mAG, mBR, mPi), anova)
#------------------------------------------------------------------------------------------
# agora faremos os gráficos dos valores observados com as retas dentro
# jitter() é para "agitar" os valores de dose antes de plotar, apenas para melhorar aspecto
layout(1)
plot(indice~jitter(dose), data=sorgo, col=sorgo$cultivar)
with(subset(sorgo, cultivar==levels(sorgo$cultivar)[1]), lines(predict(mAG)~dose, col=1))
with(subset(sorgo, cultivar==levels(sorgo$cultivar)[2]), lines(predict(mBR)~dose, col=2))
with(subset(sorgo, cultivar==levels(sorgo$cultivar)[3]), lines(predict(mPi)~dose, col=3))
#------------------------------------------------------------------------------------------
Pode-se usar as funções sapply() e/ou lapply() para tornar o código mais eficiente/enxuto.
Análise intrabloco e interbloco de experimento em blocos incompletos tipo III
Para análise intrabloco pode-se usar a função agricolae::BIB.test().
palavras-chave: #bloco_incompleto, #intrabloco, #interbloco, #efeito_aleatório, #Tukey, #média_ajustada.
#------------------------------------------------------------------------------------------
# por Walmes
#------------------------------------------------------------------------------------------
# dados de um experimento em blocos incompletos tipo III (Pimentel Gomes, p.179)
bib3 <- c(1, 1, 35, 1, 2, 28, 1, 3, 27, 2, 1, 30, 2, 2, 20, 2, 4, 22,
3, 1, 28, 3, 2, 16, 3, 5, 18, 4, 1, 36, 4, 3, 29, 4, 4, 30,
5, 1, 29, 5, 3, 19, 5, 5, 22, 6, 1, 25, 6, 4, 16, 6, 5, 19,
7, 2, 26, 7, 3, 30, 7, 4, 28, 8, 2, 27, 8, 3, 29, 8, 5, 27,
9, 2, 29, 9, 4, 29, 9, 5, 27, 10, 3, 27, 10, 4, 26, 10, 5, 29)
bib3 <- as.data.frame(matrix(bib3, ncol=3, byrow=TRUE))
names(bib3) <- c("bloc","trat","resp")
bib3 <- transform(bib3, bloc=factor(bloc), trat=factor(trat))
str(bib3)
#------------------------------------------------------------------------------------------
# número de ocorrência de tratamento com blocos e número de níveis dos fatores
with(bib3, table(trat, bloc))
levels(bib3$trat)
levels(bib3$bloc)
#------------------------------------------------------------------------------------------
# análise intrabloco
#------------------------------------------------------------------------------------------
# soma de quadrados sequencial com tratamentos ajustados aos blocos (intrabloco)
m0 <- lm(resp~bloc+trat, data=bib3)
anova(m0)
#------------------------------------------------------------------------------------------
# checagem das pressuposições
par(mfrow=c(2,2))
plot(m0)
layout(1)
#------------------------------------------------------------------------------------------
# soma de quadrados sequencial com blocos ajustados aos tratamentos
m1 <- lm(resp~trat+bloc, data=bib3)
anova(m1)
#------------------------------------------------------------------------------------------
# soma de quadrados marginal, ambos fatores ajustados um ao outro
drop1(m1, scope=.~., test="F")
car::Anova(m1, type="III")
#------------------------------------------------------------------------------------------
# obtenção das médias ajustadas
# uso da opção de contraste tipo soma zero torna mais simples as operações matriciais
# nesse caso (intercept) representa a média geral
m0 <- lm(resp~bloc+trat, data=bib3,
contrast=list(bloc=contr.sum, trat=contr.sum))
#------------------------------------------------------------------------------------------
# vamos usar a matriz de contraste dos efeitos dos tratamentos e posições dos efeitos no
# vetor de estimativas
m0$contrast$trat
m0$assign
#------------------------------------------------------------------------------------------
# obtenção das médias ajustadas é o produto matricial abaixo somado à média geral
maju <- c(m0$contrast$trat%*%coef(m0)[m0$assign==2]+coef(m0)[1])
maju
#------------------------------------------------------------------------------------------
# vetor de coeficientes do contraste tr1 vs tr2 e sua estimativa
contr <- rep(c(0,0,1,-1,0), c(1, nlevels(bib3$bloc)-1, 1, 1, nlevels(bib3$trat)-3)); contr
sum(contr*coef(m0))
#------------------------------------------------------------------------------------------
# variância do contraste, como tem equilibrio, todos os contrastes tem mesma variância
v.dif <- contr%*%vcov(m0)%*%contr; v.dif
#------------------------------------------------------------------------------------------
# diferença mínima significativa pelo teste de Tukey
delta <- qtukey(0.95, nlevels(bib3$trat), df=df.residual(m0))*sqrt(0.5*v.dif); delta
#------------------------------------------------------------------------------------------
# aplicação do teste de Tukey
taju <- agricolae::order.stat(levels(bib3$trat), maju, delta)
#------------------------------------------------------------------------------------------
# gráfico de barras com as médias
taju2 <- taju$means
names(taju2) <- taju$trt
bp <- barplot(taju2,
xlab="Tratamentos", ylab="Variável resposta",
ylim=c(0, 1.2*max(taju2)))
text(bp, taju2,
labels=paste(format(taju$means,dig=3), gsub(" ","",taju$M), sep="\n"), pos=3)
mtext(3, line=2,
text="Comparação múltipla de médias", cex=1.5)
mtext(3, line=0.5,
text="Médias seguidas de mesma letra não diferem entre si pelo teste de Tukey (5%)")
box()
#------------------------------------------------------------------------------------------
# análise interbloco
#------------------------------------------------------------------------------------------
# ajuste do modelo, anova() não faz anova e sim teste de Wald sequencial para efeitos fixos
require(nlme)
mm0 <- lme(resp~trat, random=~1|bloc, data=bib3)
anova(mm0)
#------------------------------------------------------------------------------------------
# checagem da normalidade efeitos/erros aleatórios
par(mfrow=c(1,2))
qqnorm(residuals(mm0)); qqline(residuals(mm0))
qqnorm(unlist(ranef(mm0))); qqline(unlist(ranef(mm0)))
layout(1)
#------------------------------------------------------------------------------------------
# estimativa dos componetes de variância
VarCorr(mm0)
#------------------------------------------------------------------------------------------
# matriz de contraste dos tratamentos
contr <- mm0$contrasts$trat
idx <- grep("trat", names(fixef(mm0)))
#------------------------------------------------------------------------------------------
# médias ajustadas dos tratamentos
maju <- c(contr%*%fixef(mm0)[idx]+fixef(mm0)[1])
maju
#------------------------------------------------------------------------------------------
# vetor de coeficientes do contraste tr1 vs tr2 e sua estimativa
contr <- rep(c(0, 1, -1, 0), c(1, 1, 1, nlevels(bib3$trat)-3)); contr
contr%*%maju
#------------------------------------------------------------------------------------------
# variância do contraste, como tem equilibrio, todos os contrastes tem mesma variância
v.dif <- contr%*%vcov(mm0)%*%contr; v.dif
#------------------------------------------------------------------------------------------
# diferença mínima significativa pelo teste de Tukey
delta <- c(qtukey(0.95, nlevels(bib3$trat), anova(mm0)["trat","denDF"])*sqrt(0.5*v.dif))
delta
#------------------------------------------------------------------------------------------
# aplicação do teste de Tukey
agricolae::order.stat(levels(bib3$trat), maju, delta)
#------------------------------------------------------------------------------------------
Análise intrabloco e interbloco de experimento em blocos incompletos tipo II
Temporariamente sem descrição.
palavras-chave: #bloco_incompleto, #intrabloco, #interbloco, #efeito_aleatório, #Tukey, #média_ajustada.
#------------------------------------------------------------------------------------------
# por Walmes
#------------------------------------------------------------------------------------------
# dados de um experimento em blocos incompletos tipo II (Pimentel Gomes p.188)
bib2 <- c(35,1, 28,2, 32,2, 37,3, 35,3, 25,4, 28,4, 27,5, 30,5, 32,6, 24,6, 26,7,
31,7, 27,1, 38,1, 40,3, 36,3, 27,5, 23,5, 30,7, 28,7, 25,2, 26,2, 28,4,
23,4, 24,6, 28,6, 33,1, 30,1, 22,4, 27,4, 34,7, 32,7, 39,3, 33,3, 24,6,
28,6, 34,2, 29,2, 26,5, 23,5, 33,1)
bib2 <- matrix(bib2, ncol=2, byrow=TRUE)
bib2 <- as.data.frame(bib2)
bib2$grup <- gl(3,14)
bib2$bloc <- gl(7,2)
names(bib2)[1:2] <- c("resp","trat")
bib2 <- transform(bib2, trat=factor(trat), resp=resp/10)
str(bib2)
#------------------------------------------------------------------------------------------
# análise intrabloco
#------------------------------------------------------------------------------------------
# ajuste do modelo, anova sequencial e marginal
m0 <- lm(terms(resp~grup/bloc+trat, keep.order=TRUE), data=bib2,
contrast=list(grup=contr.sum, bloc=contr.sum, trat=contr.sum))
anova(m0)
car::Anova(m0, type="III")
drop1(m0, scope=.~., test="F")
#------------------------------------------------------------------------------------------
# médias ajustadas
Xtrat <- m0$contrast$trat
assi <- m0$assign
maju <- c(Xtrat%*%coef(m0)[assi==3]+coef(m0)[1])
maju
#------------------------------------------------------------------------------------------
# vetor de coeficientes do contraste tr1 vs tr2 e sua estimativa
contr <- rep(c(0,0,0,1,-1,0), c(1,sum(assi==1),sum(assi==2),1,1,sum(assi==3)-2)); contr
sum(contr*coef(m0))
#------------------------------------------------------------------------------------------
# variância do contraste, como tem equilíbrio, todos os contrastes tem mesma variância
v.dif <- contr%*%vcov(m0)%*%contr; v.dif
#------------------------------------------------------------------------------------------
# diferença mínima significativa pelo teste de Tukey
delta <- qtukey(0.95, nlevels(bib2$trat), df=df.residual(m0))*sqrt(0.5*v.dif); delta
#------------------------------------------------------------------------------------------
# aplicação do teste de Tukey
agricolae::order.stat(levels(bib2$trat), maju, delta)
#------------------------------------------------------------------------------------------
# análise interbloco
#------------------------------------------------------------------------------------------
# criar o fator bloco dentro de grupo para associar efeito aleatório
bib2$bloc.grup <- factor(paste(bib2$grup, bib2$bloc, sep="-"))
str(bib2)
#------------------------------------------------------------------------------------------
# ajuste do modelo e teste de Wald sequencial para os efeitos fixos
require(nlme)
mm0 <- lme(resp~grup+trat, random=~1|bloc.grup, data=bib2,
contrast=list(grup=contr.sum, trat=contr.sum))
anova(mm0)
#------------------------------------------------------------------------------------------
# estimativas dos componentes de variância
VarCorr(mm0)
#------------------------------------------------------------------------------------------
# médias ajustadas
Xtrat <- mm0$contrast$trat
assi <- lapply(list(grup="grup",trat="trat"), function(x){ grep(x, names(fixef(mm0))) })
maju <- c(Xtrat%*%fixef(mm0)[assi$trat]+fixef(mm0)[1])
maju
#------------------------------------------------------------------------------------------
# vetor de coeficientes do contraste tr1 vs tr2 e sua estimativa
contr <- rep(c(0,0,1,-1,0), c(1,length(assi$grup),1,1,length(assi$trat)-2)); contr
sum(contr*fixef(mm0))
#------------------------------------------------------------------------------------------
# variância do contraste, como tem equilíbrio, todos os contrastes tem mesma variância
v.dif <- contr%*%vcov(mm0)%*%contr; v.dif
#------------------------------------------------------------------------------------------
# diferença mínima significativa pelo teste de Tukey
delta <- qtukey(0.95, nlevels(bib2$trat), df=anova(mm0)["trat","denDF"])*sqrt(0.5*v.dif)
#------------------------------------------------------------------------------------------
# aplicação do teste de Tukey
agricolae::order.stat(levels(bib2$trat), maju, delta)
#------------------------------------------------------------------------------------------
Análise intrabloco e interbloco de experimento em blocos incompletos tipo I
Temporariamente sem descrição.
palavras-chave: #bloco_incompleto, #intrabloco, #interbloco, #efeito_aleatório, #Tukey, #média_ajustada.
#------------------------------------------------------------------------------------------
# dados de um experimento em blocos incompletos tipo I (Pimentel Gomes p.185)
bib1 <- c(20,1,18,2,15,3,16,4,14,5,15,6,16,7,18,8,
24,1,18,3,25,2,19,8,13,4,16,5,12,6,16,7,
23,1,17,4,26,2,18,7,15,3,17,6,13,5,16,8,
21,1,13,5,23,2,16,3,10,4,12,7,13,6,11,8,
28,1,14,6,27,2,18,4,18,3,15,8,16,5,17,7,
22,1,17,7,24,2,16,6,18,3,14,5,15,4,17,8,
23,1,15,8,21,2,13,5,15,3,12,7,13,4,16,6)
bib1 <- matrix(bib1, ncol=2, byrow=TRUE)
bib1 <- as.data.frame(bib1)
bib1 <- cbind(rept=gl(7,8), bib1)
bib1$bloc <- gl(4,2)
names(bib1) <- c("rept","resp","trat","bloc")
bib1$trat <- factor(bib1$trat)
str(bib1)
#------------------------------------------------------------------------------------------
# análise intrabloco
#------------------------------------------------------------------------------------------
# ajuste do modelo, anova sequencial e marginal, só o teste F para trat é válido
m0 <- lm(terms(resp~rept/bloc+trat, keep.order=TRUE), data=bib1,
contrast=list(rept=contr.sum, bloc=contr.sum, trat=contr.sum))
anova(m0)
car::Anova(m0, type="III")
drop1(m0, scope=.~., test="F")
#------------------------------------------------------------------------------------------
# médias ajustadas
Xtrat <- m0$contrast$trat
assi <- m0$assign
maju <- c(Xtrat%*%coef(m0)[assi==3]+coef(m0)[1])
maju
#------------------------------------------------------------------------------------------
# vetor de coeficientes do contraste tr1 vs tr2 e sua estimativa
contr <- rep(c(0,0,0,1,-1,0), c(1,sum(assi==1),sum(assi==2),1,1,sum(assi==3)-2)); contr
sum(contr*coef(m0))
#------------------------------------------------------------------------------------------
# variância do contraste, como tem equilíbrio, todos os contrastes tem mesma variância
v.dif <- contr%*%vcov(m0)%*%contr; v.dif
#------------------------------------------------------------------------------------------
# diferença mínima significativa pelo teste de Tukey
delta <- qtukey(0.95, nlevels(bib1$trat), df=df.residual(m0))*sqrt(0.5*v.dif); delta
#------------------------------------------------------------------------------------------
# aplicação do teste de Tukey
agricolae::order.stat(levels(bib1$trat), maju, delta)
#------------------------------------------------------------------------------------------
# análise interbloco
#------------------------------------------------------------------------------------------
# criar o fator bloco dentro de repetição para associar efeito aleatório
bib1$bloc.rept <- factor(paste(bib1$rept, bib1$bloc, sep="-"))
str(bib1)
#------------------------------------------------------------------------------------------
# ajuste do modelo e teste de Wald sequencial para os efeitos fixos
mm0 <- lme(resp~rept+trat, random=~1|bloc.rept, data=bib1,
contrast=list(rept=contr.sum, trat=contr.sum))
anova(mm0)
#------------------------------------------------------------------------------------------
# estimativas dos componentes de variância
VarCorr(mm0)
#------------------------------------------------------------------------------------------
# médias ajustadas
Xtrat <- mm0$contrast$trat
assi <- lapply(list(rept="rept",trat="trat"), function(x){ grep(x, names(fixef(mm0))) })
maju <- c(Xtrat%*%fixef(mm0)[assi$trat]+fixef(mm0)[1])
maju
#------------------------------------------------------------------------------------------
# vetor de coeficientes do contraste tr1 vs tr2 e sua estimativa
contr <- rep(c(0,0,1,-1,0), c(1,length(assi$rept),1,1,length(assi$trat)-2)); contr
sum(contr*fixef(mm0))
#------------------------------------------------------------------------------------------
# variância do contraste, como tem equilíbrio, todos os contrastes tem mesma variância
v.dif <- contr%*%vcov(mm0)%*%contr; v.dif
#------------------------------------------------------------------------------------------
# diferença mínima significativa pelo teste de Tukey
delta <- qtukey(0.95, nlevels(bib1$trat), df=anova(mm0)["trat","denDF"])*sqrt(0.5*v.dif)
#------------------------------------------------------------------------------------------
# aplicação do teste de Tukey
agricolae::order.stat(levels(bib1$trat), maju, delta)
#------------------------------------------------------------------------------------------
