[R-br] Desdobrar interação dupla e tripla com glm

walmes . walmeszeviani em gmail.com
Sábado Novembro 1 17:23:16 BRST 2014


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.​
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20141101/75f677d0/attachment.html>


Mais detalhes sobre a lista de discussão R-br