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