Não foi possível enviar o arquivo. Será algum problema com as permissões?
Essa é uma revisão anterior do documento!
Regressão na análise de variância
#------------------------------------------------------------------------------------------ # 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 })}) # #------------------------------------------------------------------------------------------