Importação dos dados
##-----------------------------------------------------------------------------
## Dados de experimento feito para avaliar o efeito de período de
## armazenamento, temperatura de armazenamento e embalagem de
## armazenamento na germinação de sementes de Baccharis tridentata. O
## experimento consiste de um fatorial 3 períodos (2, 4, 6 meses) x 2
## temperaturas (5, 18°C) x 2 embalagens (papel kraft, polietileno). Uma
## testemunha foi adicionada para representar à condição inicial das
## sementes no instante 0 meses. Trata-se de um fatorial 2x2x2+1. Foram
## consideradas 4 repetições por cela experimental. As respostas medidas
## foram o número de sementes germinadas em 50 e o índice de velocidade
## de germinação.
## http://www.leg.ufpr.br/~walmes/data/ivg_milleflora.txt
da <- read.table("http://www.leg.ufpr.br/~walmes/data/ivg_tridentata.txt",
header=TRUE, sep="\t")
## Conversão para fator.
da <- transform(da,
armaz=factor(armaz),
temper=factor(temper, levels=rev(levels(temper))))
str(da)
## 'data.frame': 52 obs. of 6 variables:
## $ armaz : Factor w/ 4 levels "0","2","4","6": 1 1 1 1 2 2 2 2 2 2 ...
## $ temper : Factor w/ 3 levels "NULA","5","18": 1 1 1 1 2 2 2 2 2 2 ...
## $ embal : Factor w/ 3 levels "NULA","Papel Kraft",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ rept : int 1 2 3 4 1 2 3 4 1 2 ...
## $ ngerm50: int 31 30 29 32 31 28 31 24 21 28 ...
## $ ivg : num 2.86 2.52 2.67 2.85 2.27 ...
##-----------------------------------------------------------------------------
## Delineamento.
ftable(xtabs(~temper+embal+armaz, data=da))
## armaz 0 2 4 6
## temper embal
## NULA NULA 4 0 0 0
## Papel Kraft 0 0 0 0
## Polietileno 0 0 0 0
## 5 NULA 0 0 0 0
## Papel Kraft 0 4 4 4
## Polietileno 0 4 4 4
## 18 NULA 0 0 0 0
## Papel Kraft 0 4 4 4
## Polietileno 0 4 4 4
## Renomear o nível NULO para o primeiro nível de cada um dos
## fatores. Não tem importância pois não irá comprometer a
## estimabilidade. Isso é só para facilitar a análise porque no período
## zero de aramzenamento não houve exposição à nenhuma temperatura nem
## conservação em nenhum tipo de embalagem.
levels(da$temper)[1] <- levels(da$temper)[2]
levels(da$embal)[1] <- levels(da$embal)[2]
## Os 0's indicam celas que não existem. O 4s indicam celas presentes. É
## um fatorial incompleto.
ftable(xtabs(~temper+embal+armaz, data=da))
## armaz 0 2 4 6
## temper embal
## 5 Papel Kraft 4 4 4 4
## Polietileno 0 4 4 4
## 18 Papel Kraft 0 4 4 4
## Polietileno 0 4 4 4
##-----------------------------------------------------------------------------
## Ver.
xyplot(ngerm50~armaz|embal, groups=temper, data=da,
jitter.x=TRUE, type=c("p","a"))
##-----------------------------------------------------------------------------
## Especificação e ajuste do modelo.
## Modelo considerando distribuição normal.
m0 <- lm(ngerm50/50~armaz*temper*embal, data=da)
par(mfrow=c(2,2)); plot(m0); layout(1)
## Pode-se verificar que mesmo o dado sendo de natureza discreta os
## pressupostos de normalidade não apresentam afastamento. Nem sempre
## isso pode vir a ocorrer. As probababilidades de germinação estiveram
## entre 0.5 e 0.7 (no centro), uma faixa estreita de variação de 0.2, o
## que colabora para o não afastamento dos pressupostos.
##-----------------------------------------------------------------------------
## Modelo considerando distribuição binomial.
m0 <- glm(cbind(yes=ngerm50, no=50-ngerm50)~armaz*temper*embal,
data=da, family="quasibinomial")
## Diagnóstico.
par(mfrow=c(2,2)); plot(m0); layout(1)
## Estimativa do parâmetro de dispersão.
summary(m0)$dispersion
## [1] 1.573133
## Quadro de testes baseados na deviance para termos do modelo.
anova(m0, test="F")
## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: cbind(yes = ngerm50, no = 50 - ngerm50)
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 51 96.732
## armaz 3 0.9844 48 95.748 0.2086 0.88984
## temper 1 0.5597 47 95.188 0.3558 0.55430
## embal 1 7.5308 46 87.657 4.7871 0.03473 *
## armaz:temper 2 15.9279 44 71.730 5.0625 0.01110 *
## armaz:embal 2 5.8306 42 65.899 1.8532 0.17027
## temper:embal 1 1.2140 41 64.685 0.7717 0.38508
## armaz:temper:embal 2 2.8063 39 61.879 0.8919 0.41806
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Pode-se considerar um modelo reduzido pelo abandono de termos
## considerados nulos.
m1 <- update(m0, .~armaz*temper+embal)
## Testa se o termos abandoandos tem efeito nulo.
anova(m1, m0, test="F")
## Analysis of Deviance Table
##
## Model 1: cbind(yes = ngerm50, no = 50 - ngerm50) ~ armaz + temper + embal +
## armaz:temper
## Model 2: cbind(yes = ngerm50, no = 50 - ngerm50) ~ armaz * temper * embal
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 44 71.730
## 2 39 61.879 5 9.8509 1.2524 0.3038
## Quadro de deviance.
anova(m1, test="F")
## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: cbind(yes = ngerm50, no = 50 - ngerm50)
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 51 96.732
## armaz 3 0.9844 48 95.748 0.2046 0.89269
## temper 1 0.5597 47 95.188 0.3489 0.55775
## embal 1 7.5308 46 87.657 4.6947 0.03571 *
## armaz:temper 2 15.9279 44 71.730 4.9647 0.01137 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Mesmo o modelo reduzído contém parâmetros não estimáveis porque
## tem-se um fatorial incompleto.
coef(m1)
## (Intercept) armaz2 armaz4 armaz6 temper18
## 0.44731222 -0.09042130 -0.07003754 0.21233197 -0.40401657
## embalPolietileno armaz2:temper18 armaz4:temper18 armaz6:temper18
## -0.22987084 0.66385383 0.74072689 NA
##-----------------------------------------------------------------------------
## Ajustar modelo a partir de matriz com colunas correspondentes à
## parâmetros estimáveis.
Xmodel <- model.matrix(m1)
Xmodel <- Xmodel[,!is.na(coef(m1))]
dim(Xmodel)
## [1] 52 8
ef <- c("armaz","temper","embal")
Xm <- LSmatrix(m1, effect=ef)
grid <- equallevels(attr(Xm, "grid"), da)
grid$i <- 1:nrow(grid)
## Mantém só celas presentes.
grid <- merge(x=grid,
y=unique(subset(da, select=ef)),
all.y=TRUE, sort=FALSE)
## Matriz com funções lineares para estimar às médias.
Xm <- Xm[grid$i, !is.na(coef(m1))]
dim(Xm)
## [1] 13 8
m2 <- update(m1, .~0+Xmodel)
anova(m2, m1) ## Prova de que são o mesmo modelo.
## Analysis of Deviance Table
##
## Model 1: cbind(yes = ngerm50, no = 50 - ngerm50) ~ Xmodel - 1
## Model 2: cbind(yes = ngerm50, no = 50 - ngerm50) ~ armaz + temper + embal +
## armaz:temper
## Resid. Df Resid. Dev Df Deviance
## 1 44 71.73
## 2 44 71.73 0 0
##-----------------------------------------------------------------------------
## Estimativa das médias com intervalo de confiança com cobertura
## individual.
ci <- confint(glht(m2, linfct=Xm), calpha=univariate_calpha())
grid <- cbind(grid, ci$confint)
str(grid)
## 'data.frame': 13 obs. of 7 variables:
## $ armaz : Factor w/ 4 levels "0","2","4","6": 1 2 3 4 2 3 4 2 3 4 ...
## $ temper : Factor w/ 2 levels "5","18": 1 1 1 1 2 2 2 1 1 1 ...
## $ embal : Factor w/ 2 levels "Papel Kraft",..: 1 1 1 1 1 1 1 2 2 2 ...
## $ i : int 1 2 3 4 6 7 8 10 11 12 ...
## $ Estimate: num 0.447 0.357 0.377 0.66 0.617 ...
## $ lwr : num 0.0874 0.0853 0.1053 0.3805 0.339 ...
## $ upr : num 0.807 0.629 0.649 0.939 0.894 ...
## Estimativas na escala do preditor linear.
segplot(armaz:temper:embal~lwr+upr,
centers=Estimate, data=grid, draw=FALSE,
ylab="Cela experimental", xlab=expression(eta))
## Passando os trêm pontos pelo inverso da função de ligação.
ci <- sapply(grid[,c("Estimate","lwr","upr")], m2$family$linkinv)
colnames(ci) <- paste0("p", colnames(ci))
grid <- cbind(grid, ci)
## Estimativas na escala da resposta = probabilidade.
segplot(armaz:temper:embal~plwr+pupr,
centers=pEstimate, data=grid, draw=FALSE,
ylab="Cela experimental", xlab="Probabilidade de germinação")
## xyplot(armaz:temper:embal~ngerm50/50, data=da, pch=2,
## ylab="Cela experimental", xlab="Probabilidade de germinação",
## panel=function(y, ...){
## panel.xyplot(as.integer(y)-0.1, ...)
## })+
## as.layer(
## segplot(armaz:temper:embal~plwr+pupr,
## centers=pEstimate, data=grid, draw=FALSE))
## Valores observados em triângulo.
##-----------------------------------------------------------------------------
## Não há necessidade de se testar as médias dos tipos de embalagem
## usada uma vez que o quadro de deviance apontou haver diferença entre
## níveis. Pelo gráfico, verifica-se que a germinação no Papel Kraft é
## maior do que no Polietileno. Essa diferença é a mesma, ne escala do
## preditor linear, seja qual for a temperatura ou período de
## armazenamento uma vez que esses fatores não interagem. Para estimar o
## tamanho dessa diferença pode-se fazer.
## Usando a LSmatrix, armaz tem peso 1/4, incorreto.
LSmatrix(m1, effect="embal")
## (Intercept) armaz2 armaz4 armaz6 temper18 embalPolietileno armaz2:temper18
## [1,] 1 0.25 0.25 0.25 0.5 0 0.125
## [2,] 1 0.25 0.25 0.25 0.5 1 0.125
## armaz4:temper18 armaz6:temper18
## [1,] 0.125 0.125
## [2,] 0.125 0.125
## Obtendo sem funções prontas, peso é 1/3, o correto.
a <- grid$armaz!=0
Xe <- do.call(rbind, by(Xm[a,], INDICES=grid$embal[a], FUN=colMeans))
Xe
## (Intercept) armaz2 armaz4 armaz6 temper18 embalPolietileno
## Papel Kraft 1 0.3333333 0.3333333 0.3333333 0.5 0
## Polietileno 1 0.3333333 0.3333333 0.3333333 0.5 1
## armaz2:temper18 armaz4:temper18
## Papel Kraft 0.1666667 0.1666667
## Polietileno 0.1666667 0.1666667
## De qualquer forma, seja qual for o peso, no contraste isso não tem
## relevância pois os valores se anulam.
Xc <- apc(Xe)
summary(glht(m2, linfct=Xc))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glm(formula = cbind(yes = ngerm50, no = 50 - ngerm50) ~ Xmodel -
## 1, family = "quasibinomial", data = da)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## Papel Kraft-Polietileno == 0 0.2299 0.1058 2.172 0.0299 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
## ** A estimativa de diferença está na escala do preditor linear e não na
## escala da probabilidade.
##-----------------------------------------------------------------------------
## Desdobramento da interação armaz*temper.
## Matriz para fazer o desdobramento.
Xd <- by(Xm,
INDICES=with(grid, interaction(armaz, temper, drop=TRUE, sep=":")),
FUN=colMeans)
Xd <- do.call(rbind, Xd)
## Teste de hipótese de todos contra a testemunha.
Xc <- sweep(Xd[-1,], MARGIN=2, STAT=Xd[1,], FUN="-")
rownames(Xc) <- paste0(rownames(Xc), " vs Testemunha")
g0 <- summary(glht(m2, linfct=Xc), test=adjusted(type="none"))
g0
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glm(formula = cbind(yes = ngerm50, no = 50 - ngerm50) ~ Xmodel -
## 1, family = "quasibinomial", data = da)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 2:5 vs Testemunha == 0 -0.20536 0.22370 -0.918 0.359
## 4:5 vs Testemunha == 0 -0.18497 0.22379 -0.827 0.409
## 6:5 vs Testemunha == 0 0.09740 0.22588 0.431 0.666
## 2:18 vs Testemunha == 0 0.05448 0.22546 0.242 0.809
## 4:18 vs Testemunha == 0 0.15174 0.22646 0.670 0.503
## 6:18 vs Testemunha == 0 -0.30662 0.22335 -1.373 0.170
## (Adjusted p values reported -- none method)
## O teste mostra que não existe diferença da testemunha para nenhuma
## combinação armaz:temper. Mas cuidado, estamos fazendo isso na média
## do efeito de embalagem e isso não faz sentido. O efeito do papel
## Kraft e do polietileno então com peso 0.5 cada e já foi visto que as
## embalagens diferem. Essa germinação obtida é uma mistura de efeitos
## que na prática não representa uma possível prática. Só teria
## significado essa média ajustada se não houvesse efeito de
## embalagem. Esse exemplo mostra como não se fazer um
## desdobramento. Agora será usado não a média de efeito das duas
## embalagens mas apenas a embalagem papel Kraft que apresentou ser
## melhor que o polietileno.
## Linhas correspondentes ao Kraft.
krf <- grid$embal=="Papel Kraft" | grid$armaz=="0"
## Matriz para fazer o desdobramento.
Xk <- Xm[krf,]
rownames(Xk) <- with(grid[krf,], paste(armaz, temper, sep=":"))
## Teste de hipótese de todos contra a testemunha.
Xc <- sweep(Xk[-1,], MARGIN=2, STAT=Xk[1,], FUN="-")
rownames(Xc) <- paste0(rownames(Xc), " vs Testemunha")
g0 <- summary(glht(m2, linfct=Xc), test=adjusted(type="none"))
g0
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glm(formula = cbind(yes = ngerm50, no = 50 - ngerm50) ~ Xmodel -
## 1, family = "quasibinomial", data = da)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 2:5 vs Testemunha == 0 -0.09042 0.23004 -0.393 0.694
## 4:5 vs Testemunha == 0 -0.07004 0.23015 -0.304 0.761
## 6:5 vs Testemunha == 0 0.21233 0.23236 0.914 0.361
## 2:18 vs Testemunha == 0 0.16942 0.23193 0.730 0.465
## 4:18 vs Testemunha == 0 0.26667 0.23296 1.145 0.252
## 6:18 vs Testemunha == 0 -0.19168 0.22964 -0.835 0.404
## (Adjusted p values reported -- none method)
## Estimativas.
grid <- grid[krf,]
ci <- confint(glht(m2, linfct=Xk), calpha=univariate_calpha())
grid <- cbind(grid, ci$confint)
grid
## armaz temper embal i Estimate lwr upr pEstimate plwr
## 1 0 5 Papel Kraft 1 0.4473122 0.08743772 0.8071867 0.6100000 0.5218455
## 2 2 5 Papel Kraft 2 0.3568909 0.08526328 0.6285186 0.5882876 0.5213029
## 3 4 5 Papel Kraft 3 0.3772747 0.10530833 0.6492410 0.5932156 0.5263028
## 4 6 5 Papel Kraft 4 0.6596442 0.38053427 0.9387541 0.6591805 0.5940020
## 5 2 18 Papel Kraft 6 0.6167282 0.33900127 0.8944551 0.6494741 0.5839479
## 6 4 18 Papel Kraft 7 0.7139850 0.43296871 0.9950013 0.6712811 0.6065823
## 7 6 18 Papel Kraft 8 0.2556276 -0.01466754 0.5259228 0.5635612 0.4963332
## pupr Estimate lwr upr
## 1 0.6915097 0.4473122 0.08743772 0.8071867
## 2 0.6521535 0.3568909 0.08526328 0.6285186
## 3 0.6568394 0.3772747 0.10530833 0.6492410
## 4 0.7188479 0.6596442 0.38053427 0.9387541
## 5 0.7098087 0.6167282 0.33900127 0.8944551
## 6 0.7300746 0.7139850 0.43296871 0.9950013
## 7 0.6285317 0.2556276 -0.01466754 0.5259228
## Passando os trêm pontos pelo inverso da função de ligação.
ci <- sapply(grid[,c("Estimate","lwr","upr")], m2$family$linkinv)
colnames(ci) <- paste0("p", colnames(ci))
grid <- cbind(grid, ci)
## Estimativas na escala da resposta = probabilidade.
segplot(armaz:temper~plwr+pupr,
centers=pEstimate, data=grid, draw=FALSE,
ylab="Cela experimental", xlab="Probabilidade de germinação")
