<div dir="ltr"><div class="gmail_default" style="font-family:trebuchet ms,sans-serif">Maurício,<br><br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif">Segue código com comentários.<br><br><span style="font-family:courier new,monospace">##-----------------------------------------------------------------------------<br>## Pacotes.<br><br>require(multcomp)<br>require(doBy)<br>require(lattice)<br>require(plyr)<br><br>##-----------------------------------------------------------------------------<br>## Leitura e reconhecimento dos dados.<br><br>dados <- read.table("dados.txt", header=TRUE, sep="\t")<br>str(dados)<br><br>ftable(xtabs(~genotipo+meio+tipo, dados))<br><br>dados$prop <- with(dados, n_ex_ox/n_explant)<br><br>xyplot(prop~tipo|genotipo, groups=meio, data=dados,<br>       type=c("p","a"))<br><br>xyplot(prop~meio|tipo, groups=genotipo, data=dados,<br>       type=c("p","a"))<br><br>## Elevada presença de celas com proporção amostral igual à zero. Isso<br>## pode ser problemático pois a proporção amostral (média da amostra) é<br>## a estimativa de p no caso de uma população sem modelo de<br>## regressão. Mas p tem que estar no (0,1), não pode ser 0 nem 1. Além<br>## do mais, o erro padrão da estimativa de p envolve p*(1-p) no<br>## denominador, então se p -> 0, seu erro padrão tende à Inf. Pode ser<br>## que a estrutura experimental permita estimar valores não na borda,<br>## mas isso não é garantido.<br><br>##-----------------------------------------------------------------------------<br>## Ajuste de modelos.<br><br>## Modelo fatorial completo (maximal).<br>m0 <- glm(cbind(good=n_ex_ox, bad=n_explant-n_ex_ox)~<br>          genotipo*meio*tipo,<br>          family="quasibinomial", data=dados)<br><br>## Diagnóstico.<br>par(mfrow=c(2,2)); plot(m0); layout(1)<br><br>anova(m0,test="F")<br>drop1(m0, test="F", scope=.~.)<br><br>## Indicativo de que a tripla pode ser abandonada.<br><br>## Modelo só com efeito de primeiro grau.<br>m1 <- update(m0, .~genotipo+meio+tipo)<br>anova(m1, m0,test="F")<br><br>## Modelo só com efeitos até segundo grau.<br>m2 <- update(m0, .~(genotipo+meio+tipo)^2)<br>anova(m2, m0, test="F")<br>anova(m2, test="F")<br><br>## Indica que a dupla genotipo:tipo pode ser abandonada.<br>m3 <- update(m0, .~meio*(genotipo+tipo))<br>anova(m3, m0, test="F")<br>anova(m3, m2, test="F")<br>anova(m3, test="F")<br><br>par(mfrow=c(2,2)); plot(m2); layout(1)<br>summary(m2)<br>summary(m2)$dispersion<br><br>##-----------------------------------------------------------------------------<br>## Comparações múltiplas.<br><br>## Usando o modelo m2.<br>M <- LSmatrix(m2, effect=c("tipo", "meio"))<br><br>## Pacote wzRfun contém essas funções por conveniência.<br>## Informação para instalação em <a href="https://github.com/walmes/wzRfun">https://github.com/walmes/wzRfun</a>.<br>require(wzRfun)<br><br>grid <- equallevels(attr(M, "grid"), dados)<br>str(grid)<br><br>## Fazer comparação de "medias" de tipo estratificado por níveis de<br>## meio.<br><br>X <- by(M, INDICES=grid$meio, FUN=as.matrix)<br>X <- lapply(X, "rownames<-", levels(grid$tipo))<br>str(X)<br><br>## lapply(X, apmc, model=m2, focus="tipo", test="none")<br>L <- lapply(X, apmc, model=m2, focus="tipo", test="fdr")<br>L<br><br>L <- ldply(L)<br>names(L)[1] <- "meio"<br>L <- equallevels(L, dados)<br>str(L)<br><br>## names(m2)<br>## names(m2$family)<br><br>##-----------------------------------------------------------------------------<br>## Problemas.<br><br>L$p <- m2$family$linkinv(L$estimate)<br>xyplot(p~tipo, groups=meio, data=L, type="b")<br>xyplot(p~meio, groups=tipo, data=L, type="b")<br><br>plot(sqrt(diag(vcov(m2))))<br>summary(m2)$coefficients<br><br>## Verifica-se muito valor t próximo de zero porque o erro-padrão é<br>## demasiado grande.<br><br>## Por que o teste de comparações múltiplas não rejeitou nenhum<br>## contraste? Primeiro, o teste F é sobre uma hipótese e as comparações<br>## múltiplas são para uma série de hipóteses. Não há justificativa para<br>## que ambos concordem sempre. Nesse caso em particular, parte do<br>## resultado é porque a medida que estimativa de p vai para 0 ou 1<br>## (borda), na escala do preditor corresponde à eta indo para -Inf e Inf<br>## e o seu erro-padrão se torna cada vez maior. Um contraste tem<br>## erro-padrão que é função da matriz de covariância das estimativas e<br>## erros-apdrões altos nela implicam em erros padrões altos para os<br>## contrastes. Esse tipo de inferência, (assintótica | baseada na matriz<br>## de covariância | na aproximação quadrática da função objetivo) é<br>## chamada de Wald. Já a inferência baseada na avaliação da função<br>## objetivo (que seria a função de log-verossimilhança caso não fosse um<br>## modelo quasi) não sofre desse problema pois testa hipóteses sem<br>## depender da matriz de covariância.<br><br>## Esse resultado mostra que o planejamento de experimentos com<br>## respostas do tipo binomial deve antecipar algumas coisas para ser bem<br>## sucedido. Se espera-se que p esteja próximo das bordas, mais<br>## repetições devem ser consideradas e também um maior número de ensaios<br>## Bernoulli. Nesse experimento foram consideradas apenas 2 repetições e<br>## um número máximo de ensaios de 11. A chance de um p amostral de 0 em<br>## 11 ensaios é bem maior do que se fossem feitos 50. Talvez, se o<br>## experimento tivesse 5 repetições com 25 explantes por unidade<br>## experimental, menor a chance de tais problemas. Sabe-se que elevar o<br>## número de unidades aumenta o custo e tempo do experimento, porém, se<br>## ele não atender requisitos mínimos, custo e tempo serão demandados e<br>## as perguntas do experimentos não serão respondidas.<br><br>## Pelos gráficos, tem-se uma sugestão de que o meio M1 difere dos<br>## demais e que estes tem pouca/nenhuma diferença entre si. Não é<br>## recomendado criar hipóteses com a análise dos dados e sim antes da<br>## coleta dos mesmos, mas será explorado fazer fusão de níveis até<br>## chegar em um modelo reduzído satisfatório já que o procedimento de<br>## comparações multiplas não foi bem sucedido.<br><br>##-----------------------------------------------------------------------------<br>## Modelos reduzidos pela fusão de níveis dos fatores.<br><br>dados$meio2 <- dados$meio<br>levels(dados$meio2)<br>levels(dados$meio2) <-<br>    c("M5:M7", "M4", "M6", "M2:M3", "M1")[c(5,4,4,2,1,3,1)]<br>xtabs(~meio+meio2, dados)<br><br>m4 <- glm(cbind(good=n_ex_ox, bad=n_explant-n_ex_ox)~<br>          genotipo*meio2*tipo,<br>          family="quasibinomial", data=dados)<br>anova(m0, m4, test="F")<br>anova(m4, test="F")<br><br>dados$tipo2 <- dados$tipo<br>levels(dados$tipo2)<br>levels(dados$tipo2) <-<br>    c("T7", "T5", "T8:T6", "T3", "T2:T4", "T1" )[c(6,5,4,5,2,3,1,3)]<br>xtabs(~tipo+tipo2, dados)<br><br>m5 <- glm(cbind(good=n_ex_ox, bad=n_explant-n_ex_ox)~<br>          genotipo*meio2*tipo2,<br>          family="quasibinomial", data=dados)<br>anova(m0, m5, test="F")<br>anova(m5, test="F")<br><br>##-----------------------------------------------------------------------------<br>## Visualizando as estimativas.<br><br>M <- LSmatrix(m5, effect=c("tipo2", "meio2", "genotipo"))<br>grid <- equallevels(attr(M, "grid"), dados)<br>str(grid)<br><br>str(M)<br>ci <- confint(glht(m5, linfct=M), calpha=univariate_calpha())<br>grid <- cbind(grid, ci$confint)<br>str(grid)<br><br>grid$p <- m5$family$linkinv(grid$Estimate)<br>xyplot(p~tipo2|genotipo, groups=meio2, data=grid, type="b")<br>xyplot(p~meio2|genotipo, groups=tipo2, data=grid, type="b")<br>xyplot(p~tipo2|meio2, groups=genotipo, data=grid, type="b")<br></span><br><br></div>​À<div class="gmail_default" style="font-family:trebuchet ms,sans-serif;display:inline">​ disposição.<br></div><div class="gmail_default" style="font-family:trebuchet ms,sans-serif;display:inline">Walmes.​</div></div>