pessoais:walmes:ridiculas

Não foi possível enviar o arquivo. Será algum problema com as permissões?

Você está aqui: start » pessoais » walmes » ridiculas

Diferenças

Aqui você vê as diferenças entre duas revisões dessa página.

--- pessoais:walmes:ridiculas [2010/12/20 20:07]
walmes
+++ pessoais:walmes:ridiculas [2011/02/11 17:01] (atual)
walmes
@@ Linha 1: / Linha 1: @@
 ====== Ridiculas - dicas curtas em R ======
-===== Ridiculas - dicas curtas em R =====
+===== Regressão na análise de variância =====
+<code R>
+#------------------------------------------------------------------------------------------
+# 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
+                })})
+                                                                                          #
+#------------------------------------------------------------------------------------------
+</code>
+===== Experimento com dois fatores de efeito aditivo e perda de muitas parcelas =====
+<code R>
+#------------------------------------------------------------------------------------------
+# 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,
+.4, NA, NA, NA, NA, 258.8, 265.6, NA, NA, NA, NA, NA, 250, NA, 278.8,
+.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 de resíduos
+par(mfrow=c(2,2))
+plot(m0)
+layout(1)
+#------------------------------------------------------------------------------------------
+# estimativas dos efeitos sob a restrição do R
+summary(m0)
+#------------------------------------------------------------------------------------------
+# obtenção das médias ajustadas para os 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 multipla de médias
+require(multcomp)
+summary(glht(m0, linfct=mcp(tr="Tukey")))
+#------------------------------------------------------------------------------------------
+</code>

Diferenças

Navegação

Busca

Ferramentas

QR Code