Ajuste de modelos lineares e mistos no ambiente R

09 e 10 de Outubro de 2014 - Piracicaba - SP
Prof. Dr. Walmes M. Zeviani
Escola Superior de Agricultura “Luiz de Queiroz” - USP
Lab. de Estatística e Geoinformação - LEG
Pós Graduação em Genética e Melhoramento de Plantas Departamento de Estatística - UFPR

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))

plot of chunk unnamed-chunk-3

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

plot of chunk unnamed-chunk-3

##-----------------------------------------------------------------------------
## 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))

plot of chunk unnamed-chunk-4

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

plot of chunk unnamed-chunk-4

## 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")