Não foi possível enviar o arquivo. Será algum problema com as permissões?
Diferenças
Aqui você vê as diferenças entre duas revisões dessa página.
| - | Próxima revisão |
| - | pessoais:walmes:ridiculas [2010/12/20 20:06] walmes criada |
|---|---|
| + | Revisão anterior |
| + | pessoais:walmes:ridiculas [2011/02/11 17:01] (atual) walmes |
| @@ Linha -1,2 +1,161 @@ | |
| ====== 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, | |
| 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 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> | |
