Análise de experimento em parcela subdividida no tempo
##-----------------------------------------------------------------------------
## Definições da sessão.
## Lista de pacotes a serem usados na sessão.
pkg <- c("lattice", "latticeExtra", "gridExtra", "doBy", "multcomp",
"reshape", "plyr", "nlme", "wzRfun")
sapply(pkg, require, character.only=TRUE)
## lattice latticeExtra gridExtra doBy multcomp reshape
## TRUE TRUE TRUE TRUE TRUE TRUE
## plyr nlme wzRfun
## TRUE TRUE TRUE
## source("http://dl.dropboxusercontent.com/u/48140237/bandas.R")
sessionInfo()
## R version 3.1.1 (2014-07-10)
## Platform: x86_64-pc-linux-gnu (64-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] splines grid stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] wzRfun_0.2 nlme_3.1-117 plyr_1.8.1 reshape_0.8.5
## [5] multcomp_1.3-3 TH.data_1.0-3 mvtnorm_0.9-9997 doBy_4.5-10
## [9] MASS_7.3-34 survival_2.37-7 gridExtra_0.9.1 latticeExtra_0.6-26
## [13] RColorBrewer_1.0-5 lattice_0.20-29 knitr_1.6
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_0.10 lme4_1.1-6 Matrix_1.1-4
## [5] methods_3.1.1 minqa_1.2.3 Rcpp_0.11.2 RcppEigen_0.3.2.0.2
## [9] sandwich_2.3-0 stringr_0.6.2 tools_3.1.1 zoo_1.7-10
trellis.device(color=FALSE)
##-----------------------------------------------------------------------------
## Leitura dos dados.
da <- read.table("http://www.leg.ufpr.br/~walmes/data/soja_planta_tcc.txt",
header=TRUE, sep="\t")
str(da)
## 'data.frame': 40 obs. of 9 variables:
## $ sistema : Factor w/ 2 levels "convencional",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ adubpk : int 40 40 40 40 60 60 60 60 90 90 ...
## $ bloco : Factor w/ 4 levels "I","II","III",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ alt15.12: num 50.6 48.3 51.6 43.3 49 50.3 51 49 48.6 56.3 ...
## $ alt05.01: num 69.6 59.6 68.3 60 63 69.3 64.3 67 58 61.6 ...
## $ alt26.01: num 76.6 74.6 70.3 68 68.3 74 67 64.6 66.3 67.6 ...
## $ alt1vag : num 12.6 14 14 15.3 12.3 13.6 15.3 12.3 15 14 ...
## $ prod : num 2044 2383 3137 2889 2308 ...
## $ p100 : num 32 30.3 33.7 29 29.4 30.9 30.3 31.3 32 30.8 ...
##-----------------------------------------------------------------------------
## Empilhar os dados de altura de planta.
## Dados empilhados.
db <- melt(subset(da, select=1:6), id.vars=1:3)
## Diferença de dias.
d <- diff(as.Date(c("2006-12-15","2007-01-05","2007-01-26"))); d
## Time differences in days
## [1] 21 21
## Substituir rótulos pelo número de dias.
d <- cumsum(c(0, d)); d
## [1] 0 21 42
levels(db$variable)
## [1] "alt15.12" "alt05.01" "alt26.01"
db$dia <- d[as.integer(db$variable)]
##-----------------------------------------------------------------------------
## Ver.
xyplot(value~dia|sistema, groups=adubpk, data=db, type=c("p","a"),
auto.key=TRUE)
##-----------------------------------------------------------------------------
## O que fazer com a série de altura (de apenas 3 pontos)?
## 1) Analisar como está, efeito quadrático ou categórico;
## 2) Converter a série em média <=> área abaixo da curva ou outra.
##-----------------------------------------------------------------------------
## Ajuste do modelo de parcela subdivididas (assume independência entre
## observações de uma mesma parcela no tempo).
db <- transform(db,
adub=factor(adubpk),
D=factor(dia),
parc=interaction(bloco, sistema),
subp=interaction(bloco, sistema, adubpk))
str(db)
## 'data.frame': 120 obs. of 10 variables:
## $ sistema : Factor w/ 2 levels "convencional",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ adubpk : int 40 40 40 40 60 60 60 60 90 90 ...
## $ bloco : Factor w/ 4 levels "I","II","III",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ variable: Factor w/ 3 levels "alt15.12","alt05.01",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ value : num 50.6 48.3 51.6 43.3 49 50.3 51 49 48.6 56.3 ...
## $ dia : num 0 0 0 0 0 0 0 0 0 0 ...
## $ adub : Factor w/ 5 levels "40","60","90",..: 1 1 1 1 2 2 2 2 3 3 ...
## $ D : Factor w/ 3 levels "0","21","42": 1 1 1 1 1 1 1 1 1 1 ...
## $ parc : Factor w/ 8 levels "I.convencional",..: 5 6 7 8 5 6 7 8 5 6 ...
## $ subp : Factor w/ 40 levels "I.convencional.40",..: 5 6 7 8 13 14 15 16 21 22 ...
m0 <- lme(value~bloco+sistema*adub*D, random=~1|parc/subp, data=db,
method="ML")
##-----------------------------------------------------------------------------
## Diagnóstico.
r <- residuals(m0, type="pearson")
f <- fitted(m0)
bp <- unlist(ranef(m0)[1])
bs <- unlist(ranef(m0)[2])
grid.arrange(nrow=3,
xyplot(r~f, type=c("p", "smooth")),
xyplot(sqrt(abs(r))~f, type=c("p", "smooth")),
qqmath(r), qqmath(bp), qqmath(bs))

## Embora esteja considerando só a parte de efeito fixo serve de
## indicação.
MASS::boxcox(lm(formula(m0), data=db))

##-----------------------------------------------------------------------------
## Modelo com a variável transformada. Refazer a análise dos resíduos.
m0 <- lme(log(value)~bloco+sistema*adub*D, random=~1|parc/subp, data=db,
method="ML")
## Teste de Wald para os termos de efeito fixo.
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 60 339618 <.0001
## bloco 3 3 4 0.1261
## sistema 1 3 7 0.0812
## adub 4 24 0 0.8381
## D 2 60 478 <.0001
## sistema:adub 4 24 1 0.4773
## sistema:D 2 60 1 0.5318
## adub:D 8 60 1 0.6276
## sistema:adub:D 8 60 1 0.5189
## Nenhuma interação com os dias e nem efeito principal dos fatores
## aplicados (sistema e adubação).
##-----------------------------------------------------------------------------
## Abandono de termos.
## Teste de razão de verossimilhanças para o efeito de sistema.
m1 <- update(m0, .~bloco+sistema+adub+D)
anova(m1, m0)
## Model df AIC BIC logLik Test L.Ratio p-value
## m1 1 14 -326.3 -287.3 177.2
## m0 2 36 -305.1 -204.7 188.5 1 vs 2 22.77 0.4149
anova(m1)
## numDF denDF F-value p-value
## (Intercept) 1 78 378330 <.0001
## bloco 3 3 5 0.1110
## sistema 1 3 7 0.0719
## adub 4 28 0 0.8103
## D 2 78 478 <.0001
summary(m1)
## Linear mixed-effects model fit by maximum likelihood
## Data: db
## AIC BIC logLik
## -326.3 -287.3 177.2
##
## Random effects:
## Formula: ~1 | parc
## (Intercept)
## StdDev: 8.565e-07
##
## Formula: ~1 | subp %in% parc
## (Intercept) Residual
## StdDev: 0.02885 0.04911
##
## Fixed effects: log(value) ~ bloco + sistema + adub + D
## Value Std.Error DF t-value p-value
## (Intercept) 3.944 0.02120 78 186.02 0.0000
## blocoII 0.035 0.01898 3 1.86 0.1601
## blocoIII 0.004 0.01898 3 0.20 0.8563
## blocoIV -0.038 0.01898 3 -1.98 0.1417
## sistemadireto -0.037 0.01342 3 -2.73 0.0719
## adub60 -0.017 0.02122 28 -0.80 0.4300
## adub90 0.005 0.02122 28 0.22 0.8302
## adub120 0.005 0.02122 28 0.25 0.8050
## adub150 0.004 0.02122 28 0.21 0.8362
## D21 0.270 0.01152 78 23.44 0.0000
## D42 0.336 0.01152 78 29.17 0.0000
## Correlation:
## (Intr) blocII blcIII blocIV sstmdr adub60 adub90 adb120 adb150 D21
## blocoII -0.448
## blocoIII -0.448 0.500
## blocoIV -0.448 0.500 0.500
## sistemadireto -0.317 0.000 0.000 0.000
## adub60 -0.500 0.000 0.000 0.000 0.000
## adub90 -0.500 0.000 0.000 0.000 0.000 0.500
## adub120 -0.500 0.000 0.000 0.000 0.000 0.500 0.500
## adub150 -0.500 0.000 0.000 0.000 0.000 0.500 0.500 0.500
## D21 -0.272 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
## D42 -0.272 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.500
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.00834 -0.60959 -0.03931 0.56037 2.37675
##
## Number of Observations: 120
## Number of Groups:
## parc subp %in% parc
## 8 40
##-----------------------------------------------------------------------------
## O detalhe da parcela subdivida no tempo (ou profundidade) é que pelo
## fato de não haver aleatorização e o tempo ser uma série, pode haver
## uma correlação entre observações dependente da distância entre elas
## ao longo da série. Ou seja, medidas mais proximas são mais
## correlacionadas. Pode-se declarar esse modelo e verificar o quanto
## essa correlação significa.
## Níveis de dia em 1, 2 e 3 (são equidistantes).
db$diai <- as.integer(db$D)
m2 <- update(m0, correlation=corAR1(form=~diai|parc/subp))
## m2 <- update(m0, correlation=corSymm(form=~diai|parc/subp))
summary(m2)
## Linear mixed-effects model fit by maximum likelihood
## Data: db
## AIC BIC logLik
## -303.8 -200.7 188.9
##
## Random effects:
## Formula: ~1 | parc
## (Intercept)
## StdDev: 8.941e-07
##
## Formula: ~1 | subp %in% parc
## (Intercept) Residual
## StdDev: 0.01937 0.0484
##
## Correlation Structure: AR(1)
## Formula: ~diai | parc/subp
## Parameter estimate(s):
## Phi
## 0.221
## Fixed effects: log(value) ~ bloco + sistema * adub * D
## Value Std.Error DF t-value p-value
## (Intercept) 3.972 0.03290 60 120.72 0.0000
## blocoII 0.037 0.01968 3 1.86 0.1603
## blocoIII 0.002 0.01968 3 0.10 0.9290
## blocoIV -0.037 0.01968 3 -1.89 0.1545
## sistemadireto -0.094 0.04329 3 -2.16 0.1190
## adub60 -0.050 0.04329 24 -1.16 0.2580
## adub90 -0.048 0.04329 24 -1.12 0.2737
## adub120 0.029 0.04329 24 0.68 0.5056
## adub150 -0.033 0.04329 24 -0.77 0.4470
## D21 0.217 0.03547 60 6.12 0.0000
## D42 0.317 0.03920 60 8.08 0.0000
## sistemadireto:adub60 0.080 0.06122 24 1.31 0.2023
## sistemadireto:adub90 0.111 0.06122 24 1.81 0.0830
## sistemadireto:adub120 -0.034 0.06122 24 -0.55 0.5853
## sistemadireto:adub150 0.056 0.06122 24 0.92 0.3671
## sistemadireto:D21 0.067 0.05017 60 1.33 0.1880
## sistemadireto:D42 0.086 0.05544 60 1.55 0.1265
## adub60:D21 0.040 0.05017 60 0.79 0.4328
## adub90:D21 0.093 0.05017 60 1.85 0.0696
## adub120:D21 0.007 0.05017 60 0.13 0.8971
## adub150:D21 0.082 0.05017 60 1.63 0.1088
## adub60:D42 0.009 0.05544 60 0.16 0.8710
## adub90:D42 0.035 0.05544 60 0.64 0.5268
## adub120:D42 -0.014 0.05544 60 -0.25 0.8064
## adub150:D42 0.003 0.05544 60 0.06 0.9531
## sistemadireto:adub60:D21 -0.044 0.07095 60 -0.63 0.5335
## sistemadireto:adub90:D21 -0.151 0.07095 60 -2.13 0.0370
## sistemadireto:adub120:D21 0.003 0.07095 60 0.04 0.9674
## sistemadireto:adub150:D21 -0.053 0.07095 60 -0.74 0.4619
## sistemadireto:adub60:D42 -0.095 0.07840 60 -1.21 0.2319
## sistemadireto:adub90:D42 -0.118 0.07840 60 -1.51 0.1367
## sistemadireto:adub120:D42 -0.031 0.07840 60 -0.39 0.6951
## sistemadireto:adub150:D42 -0.059 0.07840 60 -0.75 0.4562
## Correlation:
## (Intr) blocII blcIII blocIV sstmdr adub60 adub90 adb120 adb150
## blocoII -0.299
## blocoIII -0.299 0.500
## blocoIV -0.299 0.500 0.500
## sistemadireto -0.658 0.000 0.000 0.000
## adub60 -0.658 0.000 0.000 0.000 0.500
## adub90 -0.658 0.000 0.000 0.000 0.500 0.500
## adub120 -0.658 0.000 0.000 0.000 0.500 0.500 0.500
## adub150 -0.658 0.000 0.000 0.000 0.500 0.500 0.500 0.500
## D21 -0.539 0.000 0.000 0.000 0.410 0.410 0.410 0.410 0.410
## D42 -0.596 0.000 0.000 0.000 0.453 0.453 0.453 0.453 0.453
## sistemadireto:adub60 0.465 0.000 0.000 0.000 -0.707 -0.707 -0.354 -0.354 -0.354
## sistemadireto:adub90 0.465 0.000 0.000 0.000 -0.707 -0.354 -0.707 -0.354 -0.354
## sistemadireto:adub120 0.465 0.000 0.000 0.000 -0.707 -0.354 -0.354 -0.707 -0.354
## sistemadireto:adub150 0.465 0.000 0.000 0.000 -0.707 -0.354 -0.354 -0.354 -0.707
## sistemadireto:D21 0.381 0.000 0.000 0.000 -0.579 -0.290 -0.290 -0.290 -0.290
## sistemadireto:D42 0.421 0.000 0.000 0.000 -0.640 -0.320 -0.320 -0.320 -0.320
## adub60:D21 0.381 0.000 0.000 0.000 -0.290 -0.579 -0.290 -0.290 -0.290
## adub90:D21 0.381 0.000 0.000 0.000 -0.290 -0.290 -0.579 -0.290 -0.290
## adub120:D21 0.381 0.000 0.000 0.000 -0.290 -0.290 -0.290 -0.579 -0.290
## adub150:D21 0.381 0.000 0.000 0.000 -0.290 -0.290 -0.290 -0.290 -0.579
## adub60:D42 0.421 0.000 0.000 0.000 -0.320 -0.640 -0.320 -0.320 -0.320
## adub90:D42 0.421 0.000 0.000 0.000 -0.320 -0.320 -0.640 -0.320 -0.320
## adub120:D42 0.421 0.000 0.000 0.000 -0.320 -0.320 -0.320 -0.640 -0.320
## adub150:D42 0.421 0.000 0.000 0.000 -0.320 -0.320 -0.320 -0.320 -0.640
## sistemadireto:adub60:D21 -0.270 0.000 0.000 0.000 0.410 0.410 0.205 0.205 0.205
## sistemadireto:adub90:D21 -0.270 0.000 0.000 0.000 0.410 0.205 0.410 0.205 0.205
## sistemadireto:adub120:D21 -0.270 0.000 0.000 0.000 0.410 0.205 0.205 0.410 0.205
## sistemadireto:adub150:D21 -0.270 0.000 0.000 0.000 0.410 0.205 0.205 0.205 0.410
## sistemadireto:adub60:D42 -0.298 0.000 0.000 0.000 0.453 0.453 0.226 0.226 0.226
## sistemadireto:adub90:D42 -0.298 0.000 0.000 0.000 0.453 0.226 0.453 0.226 0.226
## sistemadireto:adub120:D42 -0.298 0.000 0.000 0.000 0.453 0.226 0.226 0.453 0.226
## sistemadireto:adub150:D42 -0.298 0.000 0.000 0.000 0.453 0.226 0.226 0.226 0.453
## D21 D42 sst:60 sst:90 ss:120 ss:150 ss:D21 ss:D42 a60:D2
## blocoII
## blocoIII
## blocoIV
## sistemadireto
## adub60
## adub90
## adub120
## adub150
## D21
## D42 0.552
## sistemadireto:adub60 -0.290 -0.320
## sistemadireto:adub90 -0.290 -0.320 0.500
## sistemadireto:adub120 -0.290 -0.320 0.500 0.500
## sistemadireto:adub150 -0.290 -0.320 0.500 0.500 0.500
## sistemadireto:D21 -0.707 -0.391 0.410 0.410 0.410 0.410
## sistemadireto:D42 -0.391 -0.707 0.453 0.453 0.453 0.453 0.552
## adub60:D21 -0.707 -0.391 0.410 0.205 0.205 0.205 0.500 0.276
## adub90:D21 -0.707 -0.391 0.205 0.410 0.205 0.205 0.500 0.276 0.500
## adub120:D21 -0.707 -0.391 0.205 0.205 0.410 0.205 0.500 0.276 0.500
## adub150:D21 -0.707 -0.391 0.205 0.205 0.205 0.410 0.500 0.276 0.500
## adub60:D42 -0.391 -0.707 0.453 0.226 0.226 0.226 0.276 0.500 0.552
## adub90:D42 -0.391 -0.707 0.226 0.453 0.226 0.226 0.276 0.500 0.276
## adub120:D42 -0.391 -0.707 0.226 0.226 0.453 0.226 0.276 0.500 0.276
## adub150:D42 -0.391 -0.707 0.226 0.226 0.226 0.453 0.276 0.500 0.276
## sistemadireto:adub60:D21 0.500 0.276 -0.579 -0.290 -0.290 -0.290 -0.707 -0.391 -0.707
## sistemadireto:adub90:D21 0.500 0.276 -0.290 -0.579 -0.290 -0.290 -0.707 -0.391 -0.354
## sistemadireto:adub120:D21 0.500 0.276 -0.290 -0.290 -0.579 -0.290 -0.707 -0.391 -0.354
## sistemadireto:adub150:D21 0.500 0.276 -0.290 -0.290 -0.290 -0.579 -0.707 -0.391 -0.354
## sistemadireto:adub60:D42 0.276 0.500 -0.640 -0.320 -0.320 -0.320 -0.391 -0.707 -0.391
## sistemadireto:adub90:D42 0.276 0.500 -0.320 -0.640 -0.320 -0.320 -0.391 -0.707 -0.195
## sistemadireto:adub120:D42 0.276 0.500 -0.320 -0.320 -0.640 -0.320 -0.391 -0.707 -0.195
## sistemadireto:adub150:D42 0.276 0.500 -0.320 -0.320 -0.320 -0.640 -0.391 -0.707 -0.195
## a90:D2 a120:D2 a150:D2 a60:D4 a90:D4 a120:D4 a150:D4 s:60:D2
## blocoII
## blocoIII
## blocoIV
## sistemadireto
## adub60
## adub90
## adub120
## adub150
## D21
## D42
## sistemadireto:adub60
## sistemadireto:adub90
## sistemadireto:adub120
## sistemadireto:adub150
## sistemadireto:D21
## sistemadireto:D42
## adub60:D21
## adub90:D21
## adub120:D21 0.500
## adub150:D21 0.500 0.500
## adub60:D42 0.276 0.276 0.276
## adub90:D42 0.552 0.276 0.276 0.500
## adub120:D42 0.276 0.552 0.276 0.500 0.500
## adub150:D42 0.276 0.276 0.552 0.500 0.500 0.500
## sistemadireto:adub60:D21 -0.354 -0.354 -0.354 -0.391 -0.195 -0.195 -0.195
## sistemadireto:adub90:D21 -0.707 -0.354 -0.354 -0.195 -0.391 -0.195 -0.195 0.500
## sistemadireto:adub120:D21 -0.354 -0.707 -0.354 -0.195 -0.195 -0.391 -0.195 0.500
## sistemadireto:adub150:D21 -0.354 -0.354 -0.707 -0.195 -0.195 -0.195 -0.391 0.500
## sistemadireto:adub60:D42 -0.195 -0.195 -0.195 -0.707 -0.354 -0.354 -0.354 0.552
## sistemadireto:adub90:D42 -0.391 -0.195 -0.195 -0.354 -0.707 -0.354 -0.354 0.276
## sistemadireto:adub120:D42 -0.195 -0.391 -0.195 -0.354 -0.354 -0.707 -0.354 0.276
## sistemadireto:adub150:D42 -0.195 -0.195 -0.391 -0.354 -0.354 -0.354 -0.707 0.276
## s:90:D2 s:120:D2 s:150:D2 s:60:D4 s:90:D4 s:120:D4
## blocoII
## blocoIII
## blocoIV
## sistemadireto
## adub60
## adub90
## adub120
## adub150
## D21
## D42
## sistemadireto:adub60
## sistemadireto:adub90
## sistemadireto:adub120
## sistemadireto:adub150
## sistemadireto:D21
## sistemadireto:D42
## adub60:D21
## adub90:D21
## adub120:D21
## adub150:D21
## adub60:D42
## adub90:D42
## adub120:D42
## adub150:D42
## sistemadireto:adub60:D21
## sistemadireto:adub90:D21
## sistemadireto:adub120:D21 0.500
## sistemadireto:adub150:D21 0.500 0.500
## sistemadireto:adub60:D42 0.276 0.276 0.276
## sistemadireto:adub90:D42 0.552 0.276 0.276 0.500
## sistemadireto:adub120:D42 0.276 0.552 0.276 0.500 0.500
## sistemadireto:adub150:D42 0.276 0.276 0.552 0.500 0.500 0.500
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.03408 -0.70012 0.05905 0.66197 3.17993
##
## Number of Observations: 120
## Number of Groups:
## parc subp %in% parc
## 8 40
## Teste para correlação serial do tipo AR1 (intervalos equidistantes).
anova(m0, m2)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 36 -305.1 -204.7 188.5
## m2 2 37 -303.8 -200.7 188.9 1 vs 2 0.743 0.3887
anova(m0)
## numDF denDF F-value p-value
## (Intercept) 1 60 339618 <.0001
## bloco 3 3 4 0.1261
## sistema 1 3 7 0.0812
## adub 4 24 0 0.8381
## D 2 60 478 <.0001
## sistema:adub 4 24 1 0.4773
## sistema:D 2 60 1 0.5318
## adub:D 8 60 1 0.6276
## sistema:adub:D 8 60 1 0.5189
anova(m2)
## numDF denDF F-value p-value
## (Intercept) 1 60 350975 <.0001
## bloco 3 3 5 0.1180
## sistema 1 3 7 0.0768
## adub 4 24 0 0.8256
## D 2 60 427 <.0001
## sistema:adub 4 24 1 0.4218
## sistema:D 2 60 1 0.5742
## adub:D 8 60 1 0.5728
## sistema:adub:D 8 60 1 0.5124
##-----------------------------------------------------------------------------
## **Cuidado!** Ás vezes uma correlação pode ocorrer por uma má
## especificação da parte de média.
## Dia como linear (sendo que o efeito é não linear).
m3 <- update(m0, .~bloco+sistema*adub*dia,
correlation=corAR1(form=~diai|parc/subp))
anova(m0, m3)
## Model df AIC BIC logLik Test L.Ratio p-value
## m0 1 36 -305.1 -204.7 188.5
## m3 2 27 -254.4 -179.2 154.2 1 vs 2 68.65 <.0001
r3 <- residuals(m3, type="pearson")
f3 <- fitted(m3)
xyplot(r~f|db$subp)+layer(panel.abline(h=0, lty=2))

xyplot(r3~f3|db$subp)+layer(panel.abline(h=0, lty=2))

## methods(class="lme")
## Auto correlação.
ACF(m0)
## lag ACF
## 1 0 1.0000
## 2 1 -0.1751
## 3 2 -0.3371
## Classes de correlação.
## help(corClasses, help_type="html")