# R code from vignette source 'cap13fatadi.Rnw' # Encoding: UTF-8 #=========================================================================================== # code chunk number 1: cap13fatadi.Rnw:3-4 #=========================================================================================== options(prompt=" ", continue=" ") #=========================================================================================== # code chunk number 2: cap13fatadi.Rnw:23-57 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # experimento em fitopatologia testanto óleos essenciais e concetrações de aplicação, # foi incluída uma testemunha que serve para comparar o efeito dos óleos, até então # desconhecido sobre a ?severidade da doença. # dados (segredo está em como montar a planilha, para provocar o confundimento correto) #fa <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/fat-adi.txt", header=TRUE) fa <- read.table("../dados/fat-adi.txt", header=TRUE) str(fa) names(fa) <- c("tr","or","cc","bl","m") # nomes curtos fa <- transform(fa, cc=factor(cc), bl=factor(bl)) str(fa) fa subset(fa, or=="T" | cc==0) # veja que T e 0 ocorrem nas mesmas linhas, estão confundidos #------------------------------------------------------------------------------------------ # análise de variância para os tratamentos que despreza estrutura fatorial-adicional, nessa # fazer vamos chegar os resíduos m0 <- lm(m~bl+tr, fa) par(mfrow=c(2,2)); plot(m0); layout(1) anova(m0) #------------------------------------------------------------------------------------------ # é comum nesses experimentos o uso de notas/escalar de medida. Deve-se tomar o cuidado de # garantir a independencia entre unidades experimentais. As notas, por mais que se evite, # são comparativas. Outro ponto seria o experimento ser conduzído "as cegas", ou seja, sem # o avaliador saber qual tratamento foi aplicado a unidade experimental avaliada. Usar a # média como resposta é importante visto que as notas são discretas (ex 1 à 5) e as médias, # a medida que se tem muitas plantas por tratamento, tende a ser contínua. Pode-se adotar # uma escala contínua fazendo um risco em um segmento não graduado em uma folha de papel. # É importante aleatorizar a ordem de avaliação das plantas a cada coleta. #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 3: cap13fatadi.Rnw:62-159 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # vamos desdobrar a SQ de tr em "efeitos dos níveis fatoriais" e "níveis fatoriais vs nível # adicional". Porque? Só porque é um contraste ortogonal. A hipótese com relação a # testemunha pode não ser essa, ou seja, de contrastar com a *média* dos níveis fatoriais #------------------------------------------------------------------------------------------ # tipos de testemunha: # * negativa: aquela que se sabe antemão que dará resultado ruim, queremos saber se os # níveis avaliados são melhores (hipótese unilateral), é o caso da testemunha sem capina # em experimentos com herbicida e da testemunha em aplicação na fitopatologia; # * posititva: aquela que possivelmente dará resultado melhor ou igual aos níveis estudados # é o caso da testemunha capinada à mão em experimentos com herbicida (não sofre o efeito # negativo da toxidade); # * padrão: é o tipo de manejo padrão, o experimento é feito para encontrar alguém superior # ao padrão, é o caso das cultivares já usadas em uma região, ou da dose já usada de um # produto; #------------------------------------------------------------------------------------------ # as matrizes de contrastes envolvidas contrasts(fa$tr) contrasts(fa$or) contrasts(fa$cc) #------------------------------------------------------------------------------------------ # vamos usar o contraste de Helmet que já tá pronto contrasts(C(factor(1:5), helmert)) # último nível vs os anteriores #------------------------------------------------------------------------------------------ # passar os níveis que represetam a testemunha para a última posição fa$cc <- factor(fa$cc, levels=c("25","50","75","0")) contrasts(C(fa$cc, helmert)) # nível 0 vs demais levels(fa$cc) # basta fazer em um dos fatores fa$cc <- factor(fa$cc, levels=c("25","50","75","0")) levels(fa$or) # testemunha já é o último nível #------------------------------------------------------------------------------------------ # anova só da parte fatorial, isso para vermos as SQ e conferir nosso desdobramento anova(m0) m1 <- aov(m~bl+or*cc, data=subset(fa, tr!="TEST")) # exlui a TEST summary(m1) #------------------------------------------------------------------------------------------ # anova com fornecimento dos contrates e "arrancando" a SQ do contraste com o adicional # da SQ do fator origem m2 <- aov(m~bl+or*cc, data=fa, contrast=list(or=contr.helmert, cc=contr.helmert)) summary(m2) # or tem 3 gl mas deveria ter 2, isso porque test tá junto dessa SQ summary(m2, expand.split=FALSE, split=list("or"=list("fatorial"=c(1:2), "adicional"=3))) # separa as SQs #------------------------------------------------------------------------------------------ # nessa análise a anova diz tudo: concentação não depende do óleo, as concetrações e óleos # não diferem entre sí, mas na média de seus efeitos, são diferentes da testemunha #------------------------------------------------------------------------------------------ # teste de média da test contra as origens na menor concentração (testemunha negativa) require(agricolae) glr <- df.residual(m2) qmr <- deviance(m2)/glr with(subset(fa, cc%in%c("0","25")), HSD.test(m, or, glr, qmr)) with(subset(fa, cc%in%c("0","50")), HSD.test(m, or, glr, qmr)) with(subset(fa, cc%in%c("0","75")), HSD.test(m, or, glr, qmr)) #------------------------------------------------------------------------------------------ # uma hipótese a ser testada é saber se a testemunha difere das cc na média dos or fa$cc <- factor(fa$cc, levels=c("0","25","50","75")) m3 <- lm(m~bl+cc, data=fa) summary(m3) # test difere das concentrações, mas estas não diferem entre sí par(mfrow=c(2,2)); plot(m3); layout(1) anova(m3, m2) # não tem falta de ajuste, posso fazer o abandono de termos #------------------------------------------------------------------------------------------ # ou testar se a test difere da média dos fatorais (só porque não teve efeito de nada) require(contrast) c0 <- contrast(m3, type="average", list(bl=levels(fa$bl), cc="0")) # média da testemunha print(c0, X=TRUE) c1 <- contrast(m3, type="average", list(bl=levels(fa$bl), cc=c("25","50","75"))) # média cc print(c1, X=TRUE) c2 <- contrast(m3, type="average", list(bl=levels(fa$bl)[1], cc="0"), list(bl=levels(fa$bl)[1], cc=c("25","50","75"))) print(c2, X=TRUE) #------------------------------------------------------------------------------------------ # o intervalo de confiança da diferença é muito útil, da noção do ganho em controle da # doença com uso de qualquer concentração (entre 25 e 75) de qualquer óleo. #------------------------------------------------------------------------------------------ # o que fazer nos próximos experimentos? Usar doses menores, tentar encontrar a dose # limitante. Se for para estudar fatores de níveis contínuos, use mais níveis, explore a # região experimental usando mais níveis para melhor entender/modelar a relação entre a # resposta e a covariável. Teste de médias para fatores contínuos é gastar dinheiro e tempo # com algo e não usufruir dos benefícios. Deve-se planejar. #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 4: cap13fatadi.Rnw:164-177 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # carrega o pacote require(ExpDes) fa.fat <- fa[fa$tr!="TEST",] # tem que separar nas porções fatorial fa.adi <- fa[fa$tr=="TEST",] # e adicional fat2.ad.rbd(fa.fat$or, fa.fat$cc, fa.fat$bl, fa.fat$m, fa.adi$m) # essa função compara a média da testemunha com a média do fatorial. Na maioria dos casos # será uma hipótese sem significado, como por exemplo no caso em que houver interação entre # os fatores. #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 5: cap13fatadi.Rnw:182-225 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # experimento com óleos essenciais aplicados sob duas formas onde observou a área sobre a # curve de progresso da doença. Existem dois níveis adicionais (test - e test +). Colunas # começadas em "a" são as notas datas em cada avaliação. Teremos que calculas a aacpd. ol <- read.table("../dados/oleos.txt", header=TRUE, sep="\t") str(ol) ol$bloco <- factor(ol$bloco) #------------------------------------------------------------------------------------------ # análise gráfica, passar dados para o formato longo require(reshape) aux <- melt(ol, id.vars=1:3, variable_name="avaliacao") str(aux) #------------------------------------------------------------------------------------------ # gráfico require(lattice) xyplot(value~avaliacao|oleo, groups=forma, data=aux, type=c("p","a")) #------------------------------------------------------------------------------------------ # obter a área abaixo da curva de progresso da doença daval <- c("2010-05-13","2010-05-20","2010-05-27","2010-06-03","2010-06-10") # datas aval daval <- as.Date(daval) t <- cumsum(c(0, diff(daval))) # dias após primeira avaliação t # foi semanal #------------------------------------------------------------------------------------------ # fazendo o calculo da audpc com a agricolae::audpc() e plyr::ddply() require(agricolae) ol$audpc <- 0 for(i in 1:nrow(ol)){ ol$audpc[i] <- audpc(evaluation=ol[i,4:8], dates=t) } str(ol) xyplot(audpc~oleo, groups=forma, data=ol, type=c("p","a")) #------------------------------------------------------------------------------------------ # ajuste do modelo m0 <- lm(audpc~bloco+oleo*forma, data=ol) par(mfrow=c(2,2)); plot(m0); layout(1) anova(m0) # se deu interação não preciso nem desdobrar para ver os efeitos principais #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 6: cap13fatadi.Rnw:230-273 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # farei o desdobramento da SQ para apresentar o procedimento # ordem dos níveis deve ser mudada, deixar assim A, C, H, fung(positiva), test(negativa) levels(ol$oleo) ol$oleo <- factor(ol$oleo, levels=c("A","C","H","fung","test")) levels(ol$oleo) levels(ol$forma) ol$forma <- factor(ol$forma, levels=c("extrato","oleo","fung","test")) levels(ol$forma) #------------------------------------------------------------------------------------------ # fazer a análise usando contr.helmert porque faz com que último confronte com os demais m1 <- aov(audpc~bloco+oleo*forma, data=ol, contrast=list(oleo=contr.helmert, forma=contr.helmert)) summary(m1, split=list(oleo=list("fat"=1:2, "fung"=3, "test"=4))) summary(m1, expand.split=FALSE, split=list(oleo=list("fat"=1:2, "fung"=3, "test"=4))) # dado que existe interação, nem olhe para as SQ de efeitos principais, não tem significado #------------------------------------------------------------------------------------------ # desdobrar a interação dentro da anova, forma dentro de oleo m2 <- aov(audpc~bloco+oleo/forma, data=ol, contrast=list(oleo=contr.helmert, forma=contr.helmert)) summary(m2) coef(m2) summary(m2, expand.split=FALSE, split=list( "oleo"=list("fung"=3, "test"=4, "fat"=1:2), "oleo:forma"=list(A=1, B=2, C=3))) #------------------------------------------------------------------------------------------ # desdobrar a interação dentro da anova, oleo dentro de forma m3 <- aov(audpc~bloco+forma/oleo, data=ol, contrast=list(oleo=contr.helmert, forma=contr.helmert)) summary(m3) coef(m3) summary(m3, expand.split=FALSE, split=list( "forma"=list("fung"=2, "test"=3, "fat"=1), "forma:oleo"=list(extrato=c(1,3), oleo=c(2,4)))) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 7: cap13fatadi.Rnw:278-313 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # aplicando testes de médias ignorando a estrutura fatorial, visto que houve interação require(agricolae) ol$tr <- with(ol, interaction(oleo, forma)) glr <- df.residual(m0) qmr <- deviance(m0)/glr with(ol, HSD.test(audpc, tr, glr, qmr)) #------------------------------------------------------------------------------------------ # fazendo os contrastes m4 <- aov(audpc~bloco+oleo*forma, data=ol, contrast=list(oleo=contr.helmert, forma=contr.helmert)) summary.lm(m4) # omite os NA e a gmodels::estimable() pode ser usada # como escrever os contrastes nessa parametrização???? require(gmodels) estimable(m4, rep(c(1,0), c(1,9))) #------------------------------------------------------------------------------------------ # mudar para uma opção mais fácil: ajustar o modelo de efeito de médias levels(ol$tr) ol$tr <- factor(ol$tr, levels=rev(levels(ol$tr))) # faz a testemunha ser o primeiro nível m5 <- lm(audpc~bloco+tr, data=ol) coef(m5) # aqui é mais fácil estabelecer os contrastes, podemos usar a contrast() #------------------------------------------------------------------------------------------ # fazendo contrastes contra a testemunha negativa e positiva tr.lev <- names(coef(m5))[m5$assign==2] tr.lev <- gsub("^tr", "", tr.lev) tr.lev contrast(m5, list(bloco="1", tr="test.test"), list(bloco="1", tr=tr.lev)) contrast(m5, list(bloco="1", tr=tr.lev[1]), list(bloco="1", tr=tr.lev[-1])) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 8: cap13fatadi.Rnw:318-374 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # dados de mistura de duas componentes: restrição que a soma dá um, predição é para x # no intervalo unitário, reparemetrização torna mais clara a interpretação mis <- read.table("../dados/biomassa.txt", header=TRUE, sep="\t") str(mis) names(mis) <- tolower(names(mis)) # passando nomes das colunas para minúsculo #------------------------------------------------------------------------------------------ # por enquanto vamos análisar apenas uma espécie mis <- subset(mis, pla=="E") str(mis) #------------------------------------------------------------------------------------------ # fazendo o gráficos de todas as respostas require(reshape) aux <- melt(mis, id.vars="k", measure.vars=8:20, variable_name="resposta") # empilhamento str(aux) #------------------------------------------------------------------------------------------ # gráfico require(lattice) xyplot(value~k|resposta, data=aux, type=c("p","a","smooth","r"), scales="free") #------------------------------------------------------------------------------------------ # análise do modelo de regressão quadrático em K com parametrização usual e adicionais str(mis) m0 <- lm(dc~bj+naplus+k+I(k^2), data=mis) summary(m0) par(mfrow=c(2,2)); plot(m0); layout(1) #------------------------------------------------------------------------------------------ # análise do modelo de regressão (quadrático) mas com a parametrização de mistura m1 <- lm(dc~-1+k+na+k:na+bj+naplus, data=mis) m1 <- lm(dc~-1+k+I(1-k)+k:I(1-k)+bj+naplus, data=mis) summary(m1) # vantagem é a interpretação dos parâmetros par(mfrow=c(2,2)); plot(m1); layout(1) #------------------------------------------------------------------------------------------ # gráfico do ajuste do modelo pred <- data.frame(k=seq(0,1,l=25), bj=0, naplus=0) pred <- rbind(pred, c(0,1,0), c(0,0,1)) pred$y <- predict(m1, newdata=pred, interval="confidence") mis$kaux <- mis$k mis$kaux[mis$bj==1] <- 1.3 mis$kaux[mis$naplus==1] <- 1.6 pred$kaux <- pred$k pred$kaux[pred$bj==1] <- 1.3 pred$kaux[pred$naplus==1] <- 1.6 plot(dc~kaux, data=mis) with(subset(pred, bj<1 & naplus<1), matlines(kaux, y, type="l", col=c(1,2,2), lty=c(1,2,2))) with(subset(pred, bj==1), arrows(kaux, y[,2], kaux, y[,3], code=3, angle=90, length=0.1)) with(subset(pred, naplus==1), arrows(kaux, y[,2], kaux, y[,3], code=3, angle=90, length=0.1)) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 9: cap13fatadi.Rnw:380-409 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # a dose de melhor mistura é diferente dos componentes puros? coef(m1) kmax <- (2*m1$coef[5]/(m1$coef[5]+m1$coef[1]-m1$coef[2]))^(-1) names(kmax) <- NULL abline(v=kmax) #------------------------------------------------------------------------------------------ # usando a gmodels::estimable(), detalhe: o vetor passado não pode ter atributo names require(gmodels) ymax <- rbind("max"=c(kmax,1-kmax,0,0,kmax*(1-kmax))) ytes <- rbind("bj"=c(0,1,1,0,0), "naplus"=c(0,1,0,1,0)) ytes-rbind(ymax, ymax) estimable(m1, cm=rbind(ymax, ytes)) # estimativas estimable(m1, cm=ytes-rbind(ymax, ymax)) # contrastes contra o máximo #------------------------------------------------------------------------------------------ # usando a contrast::contrast(), vantagem: simples de usar require(contrast) contrast(m1, list("k"=kmax, "I(1 - k)"=1-kmax, bj=0, naplus=0), list("k"=0, "I(1 - k)"=1-0, bj=0, naplus=1)) contrast(m1, list("k"=kmax, "I(1 - k)"=1-kmax, bj=0, naplus=0), list("k"=0, "I(1 - k)"=1-0, bj=1, naplus=0)) #------------------------------------------------------------------------------------------ . #=========================================================================================== # code chunk number 10: cap13fatadi.Rnw:414-431 (eval = FALSE) #=========================================================================================== . #------------------------------------------------------------------------------------------ # fazendo a casualização dos níveis dos tratamentos às unidades experimentais numeradas # usar a mesma função do DIC, já que a combinação dos níveis pode ser vista como um fator A <- gl(5, 1, labels="A") # 5 níveis | B <- gl(3, 1, labels="B") # 3 níveis | porção factorial AB <- c(outer(A, B, paste)) # 5*3 = 15 níveis | tr <- c(AB, "test-neg", "test-pos") # adiciona os níveis de testemunha r <- 6 # blocos fatadi <- design.rcbd(tr, r, kinds="Super-Duper") # se for em blocos fatadi <- design.crd(tr, r, kinds="Super-Duper") # se não for em blocos str(fatadi) write.table(fatadi, file="fatadi.txt", row.names=FALSE, sep="\t") write.table(fatadi, file="fatadi.xls", row.names=FALSE, sep="\t") #------------------------------------------------------------------------------------------ .