Experimento em arranjo fatorial duplo em DIC
##-----------------------------------------------------------------------------
## Definições da sessão, pacotes a serem usados.
pkg <- c("lattice", "latticeExtra", "doBy", "multcomp", "reshape",
"plyr", "wzRfun")
sapply(pkg, library, character.only=TRUE, logical.return=TRUE)## lattice latticeExtra doBy multcomp reshape plyr
## TRUE TRUE TRUE TRUE TRUE TRUE
## wzRfun
## TRUEsource("lattice_setup.R")
##-----------------------------------------------------------------------------
## Informações sobre as versões dos pacotes.
sessionInfo()## R version 3.1.1 (2014-07-10)
## Platform: i686-pc-linux-gnu (32-bit)
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=pt_BR.UTF-8
## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] methods splines stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.3 plyr_1.8.1 reshape_0.8.5 multcomp_1.3-7
## [5] TH.data_1.0-3 mvtnorm_0.9-9997 doBy_4.5-11 MASS_7.3-34
## [9] survival_2.37-7 latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29
## [13] rmarkdown_0.3.3 knitr_1.7
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.4 evaluate_0.5.5 formatR_1.0 grid_3.1.1 htmltools_0.2.6
## [6] Matrix_1.1-4 Rcpp_0.11.3 sandwich_2.3-0 stringr_0.6.2 tools_3.1.1
## [11] yaml_2.1.13 zoo_1.7-10## obs: Para instalar um pacote faça:
## install.packages("nome_do_pacote", dependencies=TRUE)Importação dos dados
##-----------------------------------------------------------------------------
## Lendo arquivos de dados.
## Url de um arquivo com dados.
url <- "http://www.leg.ufpr.br/~walmes/data/montgomery_14-5.txt"
## Importa os dados a partir do arquivo na web.
da <- read.table(file=url, header=TRUE, sep="\t")
names(da) <- substr(names(da), 1, 4)
## Convertendo para fator.
da <- transform(da, isol=factor(isol), temp=factor(temp))
## Mostra a estrutura do objeto.
str(da)## 'data.frame': 72 obs. of 3 variables:
## $ isol: Factor w/ 4 levels "1","2","3","4": 1 1 1 2 2 2 3 3 3 4 ...
## $ temp: Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ y : num 6.6 2.7 6 3 2.1 5.9 5.7 3.2 5.3 7 ...## Os dados são de um experimento para verificar o efeito a adubação com
## potássio nos componentes de produção de soja sob diferentes níveis de
## umidade.
##-----------------------------------------------------------------------------
## Análise exploratória.
xyplot(y~temp, groups=isol, data=da,
type=c("p","a"), auto.key=list(columns=nlevels(da$isol)))
Especificação e estimação do modelo
\[ \begin{aligned} Y|\text{fontes de variação} &\,\sim \text{Normal}(\mu_{ij}, \sigma^2) \newline \mu_{ij} &= \mu+\alpha_i+\tau_j \end{aligned} \]
em que \(\alpha_i\) é o efeito do isolante \(i\), \(\tau_j\) o efeito da temperatura \(j\) sobre a média da variável resposta \(Y\). \(\mu\) é a média de \(Y\) na ausência do efeito dos termos mencioandos. \(\mu_{ij}\) é a média para as observações na combinação \(ij\) e \(\sigma^2\) é a variância das observações ao redor desse valor médio.
##-----------------------------------------------------------------------------
## Espeficicação do modelo.
m0 <- lm(y~isol*temp, data=da)
## Resíduos.
par(mfrow=c(2,2)); plot(m0); layout(1)
## MASS::boxcox(m0)
## Quadro de anova.
anova(m0)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## isol 3 453.61 151.203 40.0658 2.421e-14 ***
## temp 2 2.44 1.222 0.3237 0.7247
## isol:temp 6 38.54 6.423 1.7019 0.1362
## Residuals 60 226.43 3.774
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##-----------------------------------------------------------------------------
## Partição da soma de quadrados.
## Essa etapa não é obrigatória porém muitos aplicativos reportam a
## particação do quadro de análise de variância como parte dos
## resultados do desdobramente da interação. Aqui será feito à titulo de
## demonstração.
## Partição para efeito de temp dentro de cada isol. Não faz muito
## sentido fazer o desdobramente porque não houve interação. Será feito
## por razões de demonstação.
m1 <- aov(y~isol/temp, data=da)
summary(m1)## Df Sum Sq Mean Sq F value Pr(>F)
## isol 3 453.6 151.20 40.066 2.42e-14 ***
## isol:temp 8 41.0 5.12 1.357 0.234
## Residuals 60 226.4 3.77
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ef <- data.frame(fv=m1$assign, par=names(coef(m1)),
stringsAsFactors=FALSE)
ef2 <- subset(ef, fv==2)
i <- paste0("^isol", levels(da$isol), ":temp")
TinI <- lapply(i, grep, x=ef2$par)
names(TinI) <- levels(da$isol)
lapply(TinI, function(i) ef2[i,])## $`1`
## fv par
## 5 2 isol1:temp2
## 9 2 isol1:temp3
##
## $`2`
## fv par
## 6 2 isol2:temp2
## 10 2 isol2:temp3
##
## $`3`
## fv par
## 7 2 isol3:temp2
## 11 2 isol3:temp3
##
## $`4`
## fv par
## 8 2 isol4:temp2
## 12 2 isol4:temp3## Quadro de análise de variância particionado.
summary(m1, split=list("isol:temp"=TinI))## Df Sum Sq Mean Sq F value Pr(>F)
## isol 3 453.6 151.20 40.066 2.42e-14 ***
## isol:temp 8 41.0 5.12 1.357 0.234
## isol:temp: 1 2 10.5 5.25 1.390 0.257
## isol:temp: 2 2 13.7 6.87 1.821 0.171
## isol:temp: 3 2 15.3 7.67 2.032 0.140
## isol:temp: 4 2 1.4 0.70 0.187 0.830
## Residuals 60 226.4 3.77
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##-----------------------------------------------------------------------------
## Outra forma de avaliar a mesma hipótese é pela estatística F do teste
## de Wald. Esse mecanismo de teste de hipótese é mais geral e pode ser
## usado em outras classes de modelo, como glm, modelos de sobrevivência
## e modelos mistos.
i <- paste0("^isol", levels(da$isol), ":temp")
TinI <- lapply(i, grep, x=ef$par)
names(TinI) <- levels(da$isol)
lapply(TinI, function(i) ef[i,])## $`1`
## fv par
## 5 2 isol1:temp2
## 9 2 isol1:temp3
##
## $`2`
## fv par
## 6 2 isol2:temp2
## 10 2 isol2:temp3
##
## $`3`
## fv par
## 7 2 isol3:temp2
## 11 2 isol3:temp3
##
## $`4`
## fv par
## 8 2 isol4:temp2
## 12 2 isol4:temp3## Matriz de covariância dos efeitos e estimativas. É um teste baseado
## em aproximação quadrática.
Sigma <- vcov(m1)
b <- coef(m1)
## Requer que o pacote 'aod' esteja instalado.
Fwald <- sapply(TinI, simplify=FALSE,
function(i){
wt <- aod::wald.test(
Sigma=Sigma,
b=b, Terms=i,
df=df.residual(m1))$result$Ftest[c(1,4)]
data.frame(as.list(wt))
})
## str(Fwald)
## Valores de F e correspondente p-valor.
do.call(rbind, Fwald)## Fstat P
## 1 1.3902649 0.2569143
## 2 1.8208584 0.1707186
## 3 2.0316652 0.1400400
## 4 0.1865169 0.8303238Comparações múltiplas
##-----------------------------------------------------------------------------
## Comparações múltiplas para isol.
## Pelo teste F não existe efeito de temperatura.
LSmeans(m0, effect="temp")## estimate se df t.stat p.value lwr upr temp
## 1 5.620833 0.3965403 60 14.17468 8.221715e-21 4.827635 6.414032 1
## 2 5.679167 0.3965403 60 14.32179 5.085369e-21 4.885968 6.472365 2
## 3 5.262500 0.3965403 60 13.27103 1.670691e-19 4.469301 6.055699 3## Para isolante, fazer comparações múltiplas.
LSmeans(m0, effect="isol")## estimate se df t.stat p.value lwr upr isol
## 1 4.183333 0.4578853 60 9.136204 5.827113e-13 3.267426 5.099240 1
## 2 2.350000 0.4578853 60 5.132290 3.252575e-06 1.434093 3.265907 2
## 3 6.511111 0.4578853 60 14.219961 7.089549e-21 5.595204 7.427018 3
## 4 9.038889 0.4578853 60 19.740509 6.284742e-28 8.122982 9.954796 4X <- LSmatrix(m0, effect="isol")
grid <- equallevels(attr(X, "grid"), da)
rownames(X) <- levels(grid$isol)
grid <- apmc(X, model=m0, focus="isol", test="single-step")
grid## isol estimate lwr upr cld
## 1 1 4.183333 3.267426 5.099240 c
## 2 2 2.350000 1.434093 3.265907 d
## 3 3 6.511111 5.595204 7.427018 b
## 4 4 9.038889 8.122982 9.954796 a##-----------------------------------------------------------------------------
## Representação gráfica dos resultados.
ylab <- "y"
xlab <- "Isolante"
sub <- list(
"Médias seguidas de mesma letra indicam diferença nula à 5%.",
font=1, cex=0.8)
## Gráfico final.
segplot(isol~lwr+upr, centers=estimate,
data=grid, draw=FALSE,
xlab=xlab, ylab=ylab, sub=sub,
letras=sprintf("%0.2f %s", grid$estimate, grid$cld),
panel=function(x, y, z, subscripts, centers, letras, ...){
panel.segplot(x=x, y=y, z=z, centers=centers,
subscripts=subscripts, ...)
panel.text(y=as.numeric(z)[subscripts],
x=centers[subscripts],
label=letras[subscripts],
srt=0, adj=c(0.5,-0.5))
})
##-----------------------------------------------------------------------------
## Análise com o pacote ExpDes.
require(ExpDes)## Loading required package: ExpDes
##
## Attaching package: 'ExpDes'
##
## The following object is masked from 'package:MASS':
##
## ginv## Funções do pacote ExpDes.
ls("package:ExpDes")## [1] "ccboot" "crd" "duncan" "fat2.ad.crd" "fat2.ad.rbd"
## [6] "fat2.crd" "fat2.rbd" "fat3.ad.crd" "fat3.ad.rbd" "fat3.crd"
## [11] "fat3.rbd" "ginv" "lastC" "latsd" "lsd"
## [16] "lsdb" "order.group" "order.stat.SNK" "rbd" "reg.poly"
## [21] "scottknott" "snk" "split2.crd" "split2.rbd" "tapply.stat"
## [26] "tukey"## Factorial (2 factors) in complete randomized design.
with(da, fat2.crd(factor1=isol, factor2=temp, resp=y,
quali=c(TRUE, TRUE), mcomp="tukey"))## ------------------------------------------------------------------------
## Legend:
## FACTOR 1: F1
## FACTOR 2: F2
## ------------------------------------------------------------------------
##
##
## Analysis of Variance Table
## ------------------------------------------------------------------------
## DF SS MS Fc Pr>Fc
## F1 3 453.61 151.203 40.066 0.00000
## F2 2 2.44 1.222 0.324 0.72471
## F1*F2 6 38.54 6.423 1.702 0.13619
## Residuals 60 226.43 3.774
## Total 71 721.02
## ------------------------------------------------------------------------
## CV = 35.19 %
##
## ------------------------------------------------------------------------
## Shapiro-Wilk normality test
## p-value: 0.4915693
## According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
## ------------------------------------------------------------------------
##
## Not significant interaction: analyzing the simple effect
## ------------------------------------------------------------------------
## F1
## Tukey's test
## ------------------------------------------------------------------------
## Groups Treatments Means
## a 4 9.038889
## b 3 6.511111
## c 1 4.183333
## d 2 2.35
## ------------------------------------------------------------------------
##
## F2
## According to the F test, the means of this factor are statistical equal.
## Levels Means
## 1 1 5.620833
## 2 2 5.679167
## 3 3 5.262500
## ------------------------------------------------------------------------## E se temperatura for considerada como numérica? Se for declarado
## efeito linear para ela?
da$tmp <- as.numeric(as.character(da$temp))
m2 <- lm(y~isol*tmp, data=da)
## Tem falta de ajuste?
anova(m2, m0)## Analysis of Variance Table
##
## Model 1: y ~ isol * tmp
## Model 2: y ~ isol * temp
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 64 228.56
## 2 60 226.43 4 2.1306 0.1411 0.9662## Quadro de anova.
anova(m2)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## isol 3 453.61 151.203 42.3385 3.375e-15 ***
## tmp 1 1.54 1.541 0.4315 0.51363
## isol:tmp 3 37.31 12.436 3.4822 0.02081 *
## Residuals 64 228.56 3.571
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Estimativas dos efeitos/coeficientes.
summary(m2)##
## Call:
## lm(formula = y ~ isol * tmp, data = da)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.0028 -1.1222 0.1014 1.1743 4.2806
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.0500 1.1785 5.134 2.87e-06 ***
## isol2 -1.5667 1.6666 -0.940 0.3507
## isol3 -1.7222 1.6666 -1.033 0.3053
## isol4 2.6056 1.6666 1.563 0.1229
## tmp -0.9333 0.5455 -1.711 0.0919 .
## isol2:tmp -0.1333 0.7715 -0.173 0.8633
## isol3:tmp 2.0250 0.7715 2.625 0.0108 *
## isol4:tmp 1.1250 0.7715 1.458 0.1497
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.89 on 64 degrees of freedom
## Multiple R-squared: 0.683, Adjusted R-squared: 0.6483
## F-statistic: 19.7 on 7 and 64 DF, p-value: 8.592e-14## Valores preditos.
da$ypred <- predict(m2)
## Obsevados vs preditos.
xyplot(y~temp, groups=isol, data=da,
auto.key=list(columns=nlevels(da$isol)))+
as.layer(
xyplot(ypred~temp, groups=isol, data=da, type="l")
)
## Quais as inclinações que diferem de zero?
m3 <- lm(y~0+isol/tmp, data=da)
summary(m3)##
## Call:
## lm(formula = y ~ 0 + isol/tmp, data = da)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.0028 -1.1222 0.1014 1.1743 4.2806
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## isol1 6.0500 1.1785 5.134 2.87e-06 ***
## isol2 4.4833 1.1785 3.804 0.000320 ***
## isol3 4.3278 1.1785 3.672 0.000492 ***
## isol4 8.6556 1.1785 7.345 4.61e-10 ***
## isol1:tmp -0.9333 0.5455 -1.711 0.091950 .
## isol2:tmp -1.0667 0.5455 -1.955 0.054921 .
## isol3:tmp 1.0917 0.5455 2.001 0.049626 *
## isol4:tmp 0.1917 0.5455 0.351 0.726489
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.89 on 64 degrees of freedom
## Multiple R-squared: 0.9216, Adjusted R-squared: 0.9118
## F-statistic: 94.05 on 8 and 64 DF, p-value: < 2.2e-16