Regressão na presença de termos de efeito aleatório

##=============================================================================
## Aplicação de modelos de regressão linear e
##   não linear em ciências agrárias
##
##   09 à 11 de Dezembro de 2014 - Goiânia/GO
##   Embrapa Arroz e Feijão
## 
##                                                  Prof. Dr. Walmes M. Zeviani
##                                                                LEG/DEST/UFPR
##=============================================================================

##-----------------------------------------------------------------------------
## Definições da sessão.

pkg <- c("lattice", "latticeExtra", "gridExtra", "car", "alr3",
         "plyr", "reshape", "doBy", "multcomp", "nlme", "rootSolve",
         "wzRfun")

sapply(pkg, require, character.only=TRUE)

##      lattice latticeExtra    gridExtra          car         alr3         plyr 
##         TRUE         TRUE         TRUE         TRUE         TRUE         TRUE 
##      reshape         doBy     multcomp         nlme    rootSolve       wzRfun 
##         TRUE         TRUE         TRUE         TRUE         TRUE         TRUE

trellis.device(color=FALSE)

Fatorial duplo analisado com regressão não linear (misto)

##-----------------------------------------------------------------------------
## Dados de experimento com soja sob 3 níveis de umidade do solo (massa
## de água/massa de solo) e níveis de potássio fornecidos na adubação. A
## hipótese sob estudo é que presença de potássio dá suporte para a
## cultura superar a ocorrência de períodos sem chuva. Experimento
## conduzido em casa de vegetação com 5 blocos, 2 plantas por
## u.e. (vaso).

## Para acessar o artigo (pdf online).
## browseURL("http://www.scielo.br/pdf/rca/v43n2/a03v43n2.pdf")

link <- "http://www.leg.ufpr.br/~walmes/data/soja.txt"
da <- read.table(link, header=TRUE, sep="\t", dec=",")
names(da) <- substr(names(da), 1, 4)

da <- subset(da[-74,], select=c("agua", "pota", "reng", "bloc"))
da$AG <- factor(da$agua)
names(da)[c(2:3)] <- c("x","y")

##-----------------------------------------------------------------------------
## Informa que bloco é de efeito aleatório.

db <- groupedData(y~x|bloc, data=da, order=FALSE)
str(db)

## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData' and 'data.frame':   74 obs. of  5 variables:
##  $ agua: num  37.5 37.5 37.5 37.5 37.5 50 50 50 50 50 ...
##  $ x   : int  0 30 60 120 180 0 30 60 120 180 ...
##  $ y   : num  14.6 21.5 24.6 21.9 28.1 ...
##  $ bloc: Factor w/ 5 levels "I","II","III",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ AG  : Factor w/ 3 levels "37.5","50","62.5": 1 1 1 1 1 2 2 2 2 2 ...
##  - attr(*, "formula")=Class 'formula' length 3 y ~ x | bloc
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##  - attr(*, "FUN")=function (x)  
##  - attr(*, "order.groups")= logi FALSE

##-----------------------------------------------------------------------------
## Visualiza os dados.

p0 <-
    xyplot(y~x|AG, data=db, layout=c(3,1), col=1,
           xlab=expression("Conteúdo de potássio no solo"~(mg~dm^{-3})),
           ylab=expression("Rendimento de grãos"~(g~vaso^{-1})))
print(p0)

##-----------------------------------------------------------------------------
## Parametrizações do modelo quadrático-platô.

require(rpanel)

## Original
##  b0: intercepto.
##  b1: coeficiente angular ou taxa na origem.
##  b2: coeficiente de concavidade/curvatura ou grau de especificidade.

fx0 <- function(panel){
    with(panel, {
        xb <- -b1/(2*b2)
        curve(b0+(b1*x+b2*x^2)*(x<=xb)+
              (b1*xb+b2*xb^2)*(x>xb), xlim[1], xlim[2], ylim=ylim)
        abline(a=b0, b=b1, lty=2)
        abline(v=xb, lty=2)
    })
    panel
}

panel <- rp.control(title="Original",
                    xlim=c(0,10), ylim=c(0,10))
rp.slider(panel, variable=b0, from=0, to=10, initval=2,
          showvalue=TRUE, action=fx0, title="b0")
rp.slider(panel, variable=b1, from=0, to=10, initval=0,
          showvalue=TRUE, action=fx0, title="b1")
rp.slider(panel, variable=b2, from=-2, to=0, initval=0,
          showvalue=TRUE, action=fx0, title="b2")

## Canônica (ver Tese Walmes).
##  yb: valor crítico em y e máximo/mínimo da função.
##  xb: valor crítico em x e ponto de troca/união.
##  b2: coeficiente de concavidade/curvatura ou grau de especificidade.

fx1 <- function(panel){
    with(panel, {
        curve(yb+(x<=xb)*b2*(x-xb)^2, xlim[1], xlim[2], ylim=ylim)
        abline(v=xb-0:1, h=yb-c(0,-b2), lty=2)
    })
    panel
}

panel <- rp.control(title="Canônica",
                    xlim=c(0,10), ylim=c(0,10))
rp.slider(panel, variable=yb, from=0, to=10, initval=5,
          showvalue=TRUE, action=fx1, title="yb")
rp.slider(panel, variable=xb, from=0, to=10, initval=5,
          showvalue=TRUE, action=fx1, title="xb")
rp.slider(panel, variable=b2, from=-5, to=0, initval=0,
          showvalue=TRUE, action=fx1, title="b2")

## Mistura original e canônica (feito para o Curso).
##  b0: intercepto.
##  b1: coeficiente angular ou taxa na origem.
##  xb: valor crítico em x e ponto de troca/união.

fx2 <- function(panel){
    with(panel, {
        curve(b0+(x<xb)*(b1*x-0.5*b1*x^2/xb)+(x>=xb)*(0.5*b1*xb),
              xlim[1], xlim[2], ylim=ylim)
        abline(v=xb-0:1, lty=2)
        abline(a=b0, b=b1, lty=2)
    })
    panel
}

panel <- rp.control(title="Mistura",
                    xlim=c(0,10), ylim=c(0,10))
rp.slider(panel, variable=b0, from=0, to=10, initval=2,
          showvalue=TRUE, action=fx2, title="b0")
rp.slider(panel, variable=b1, from=0, to=10, initval=0,
          showvalue=TRUE, action=fx2, title="b1")
rp.slider(panel, variable=xb, from=0, to=10, initval=5,
          showvalue=TRUE, action=fx2, title="xb")

##-----------------------------------------------------------------------------
## Ajuste do modelo linear-platô.

db <- groupedData(y~x|bloc, data=da, order=FALSE)

## Modelo linear platô.
n1.lps <- nlme(y~th0+th1*x*(x<=thb)+th1*thb*(x>thb),
               data=db, method="ML",
               fixed=th0+th1+thb~AG,
               random=th0~1|bloc,
               start=c(
                   15,0,0,
                   0.22,0,0,
                   40,20,50))

anova(n1.lps, type="marginal")

##                 numDF denDF   F-value p-value
## th0.(Intercept)     1    61 158.10695  <.0001
## th0.AG              2    61   1.60851  0.2086
## th1.(Intercept)     1    61  22.97344  <.0001
## th1.AG              2    61   0.07505  0.9278
## thb.(Intercept)     1    61  35.89719  <.0001
## thb.AG              2    61   8.25653  0.0007

## Modelo linear plato reduzido.
n1.lpr <- update(n1.lps,
                 fixed=list(th0+th1~1, thb~AG),
                 start=c(
                     15,
                     0.22,
                     40,20,50))

anova(n1.lps, n1.lpr)

##        Model df      AIC      BIC   logLik   Test  L.Ratio p-value
## n1.lps     1 11 347.9881 373.3328 -162.994                        
## n1.lpr     2  7 346.3440 362.4724 -166.172 1 vs 2 6.355875  0.1741

##-----------------------------------------------------------------------------
## Predição.

## summary(n1.lps)
## summary(n1.lpr)

pred <- expand.grid(x=0:180, AG=levels(da$AG))
pred$y1s <- predict(n1.lps, newdata=pred, level=0)
pred$y1r <- predict(n1.lpr, newdata=pred, level=0)

print(p0)+
    as.layer(xyplot(y1s+y1r~x|AG, data=pred, type="l"))

##-----------------------------------------------------------------------------
## Predição com bandas.

numF1 <- function(theta, xi, AG){
    theta[1]+
        (xi<=AG%*%theta[3:5])*(theta[2]*xi)+
            (xi>AG%*%theta[3:5])*(theta[2]*AG%*%theta[3:5])
}

##-----------------------------------------------------------------------------
## Faz predição com bandas.

pred1 <- expand.grid(x=seq(0,180,2), AG=levels(da$AG))
c1 <- fixef(n1.lpr)
F1 <- gradient(numF1, x=c1, xi=pred1$x,
               AG=model.matrix(~AG, pred1))
dim(F1)

## [1] 273   5

colnames(F1)==names(c1)

## [1] TRUE TRUE TRUE TRUE TRUE

head(F1)

##      th0 th1 thb.(Intercept) thb.AG50 thb.AG62.5
## [1,]   1   0               0        0          0
## [2,]   1   2               0        0          0
## [3,]   1   4               0        0          0
## [4,]   1   6               0        0          0
## [5,]   1   8               0        0          0
## [6,]   1  10               0        0          0

U1 <- chol(vcov(n1.lpr))
pred1$se <- sqrt(apply(F1%*%t(U1), 1, function(x) sum(x^2)))
zval <- qnorm(p=c(lwr=0.025, fit=0.5, upr=0.975))
me <- outer(pred1$se, zval, "*")
fit <- predict(n1.lpr, newdata=pred1, level=0)
pred1 <- cbind(pred1, sweep(me, 1, fit, "+"))
str(pred1)

## 'data.frame':    273 obs. of  6 variables:
##  $ x  : num  0 2 4 6 8 10 12 14 16 18 ...
##  $ AG : Factor w/ 3 levels "37.5","50","62.5": 1 1 1 1 1 1 1 1 1 1 ...
##  $ se : num  0.63 0.612 0.595 0.579 0.564 ...
##  $ lwr: num  13.6 14.1 14.6 15.1 15.6 ...
##  $ fit: num  14.8 15.3 15.8 16.2 16.7 ...
##  $ upr: num  16.1 16.5 16.9 17.4 17.8 ...

print(p0)+
    as.layer(xyplot(fit~x|AG, data=pred1, type="l", col=2,
                    prepanel=prepanel.cbH, cty="bands",
                    ly=pred1$lwr, uy=pred1$upr,
                    panel=panel.cbH
                    ))

##-----------------------------------------------------------------------------
## Ajusta modelo quadrático-platô reparametrizado.

## Modelo quadrático platô mistura de parametrizações.
n3.qms <- nlme(y~b0+(x<xb)*(b1*x-0.5*b1*x^2/xb)+(x>=xb)*(0.5*b1*xb),
               data=db, method="ML",
               fixed=b0+b1+xb~AG,
               random=b0~1,
               start=c(
                   15,0,0,
                   0.22,0,0,
                   88,20,40))

## n3.qms <- gnls(y~b0+(x<xb)*(b1*x-0.5*b1*x^2/xb)+(x>=xb)*(0.5*b1*xb),
##                data=db,
##                params=b0+b1+xb~AG,
##                start=c(
##                    15,0,0,
##                    0.22,0,0,
##                    88,20,40))

anova(n3.qms, type="marginal")

##                numDF denDF   F-value p-value
## b0.(Intercept)     1    61 161.33157  <.0001
## b0.AG              2    61   1.10876  0.3365
## b1.(Intercept)     1    61  24.94174  <.0001
## b1.AG              2    61   0.24038  0.7871
## xb.(Intercept)     1    61  30.87004  <.0001
## xb.AG              2    61   5.98419  0.0042

## Modelo quadrático platô mistura de parametrizações reduzido.
n3.qmr <- update(n3.qms,
                 fixed=list(b0+b1~1, xb~AG),
                 start=c(
                     15,
                     0.22,
                     88,20,40))

anova(n3.qms, n3.qmr)

##        Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## n3.qms     1 11 350.4609 375.8057 -164.2305                        
## n3.qmr     2  7 349.6555 365.7840 -167.8278 1 vs 2 7.194554   0.126

##-----------------------------------------------------------------------------
## Predição.

pred$y3s <- predict(n3.qms, newdata=pred, level=0)
pred$y3r <- predict(n3.qmr, newdata=pred, level=0)

print(p0)+
    as.layer(xyplot(y3s+y3r~x|AG, data=pred, type="l"))

##-----------------------------------------------------------------------------
## Ajusta na restrição de estimativas por casela.

n3.qmc <- update(n3.qms,
                 fixed=list(b0+b1~1, xb~AG-1),
                 start=c(
                     15,
                     0.22,
                     88,90,140))

summary(n3.qmc)$tTable

##                 Value  Std.Error DF  t-value      p-value
## b0         14.7323486  0.6748211 65 21.83149 1.245836e-31
## b1          0.2974695  0.0249935 65 11.90187 5.586195e-18
## xb.AG37.5  68.3697470  6.2813311 65 10.88460 2.746503e-16
## xb.AG50    99.6446545  8.4233734 65 11.82954 7.339089e-18
## xb.AG62.5 151.5082133 13.6860878 65 11.07024 1.337246e-16

intervals(n3.qmc)

## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                 lower        est.       upper
## b0         13.4309649  14.7323486  16.0337324
## b1          0.2492698   0.2974695   0.3456691
## xb.AG37.5  56.2562808  68.3697470  80.4832131
## xb.AG50    83.4002874  99.6446545 115.8890216
## xb.AG62.5 125.1147694 151.5082133 177.9016571
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: bloc 
##            lower      est.    upper
## sd(b0) 0.2937738 0.7802222 2.072162
## 
##  Within-group standard error:
##    lower     est.    upper 
## 1.911293 2.258342 2.668409

##-----------------------------------------------------------------------------
## Matriz gradiente, \frac{\partial \eta(x_i, \theta)}{\partial \theta_j}.

## Baseado em derivadas numéricas.
numF <- function(theta, xi, AG){
    ## b0+(x<xb)*(b1*x-0.5*b1*x^2/xb)+(x>=xb)*(0.5*b1*xb)
    theta[1]+
        (xi<AG%*%theta[3:5])*(theta[2]*xi-0.5*theta[2]*xi^2/AG%*%theta[3:5])+
            (xi>=AG%*%theta[3:5])*(0.5*theta[2]*AG%*%theta[3:5])
}

c3 <- fixef(n3.qmc)
s <- seq(0,180,30)
gradient(numF, x=c3, xi=s, AG=cbind(rep(1,length(s)),0,0))[,1:3]

##      b0       b1  xb.AG37.5
## [1,]  1  0.00000 0.00000000
## [2,]  1 23.41814 0.02863696
## [3,]  1 33.67257 0.11454785
## [4,]  1 34.18487 0.14873473
## [5,]  1 34.18487 0.14873473
## [6,]  1 34.18487 0.14873473
## [7,]  1 34.18487 0.14873473

##-----------------------------------------------------------------------------
## Faz predição com bandas.

pred3 <- expand.grid(x=seq(0,180,2), AG=levels(da$AG))

F3 <- gradient(numF, x=c3, xi=pred3$x,
               AG=model.matrix(~-1+AG, pred3))
dim(F3)

## [1] 273   5

colnames(F3)==names(c3)

## [1] TRUE TRUE TRUE TRUE TRUE

head(F3)

##      b0       b1    xb.AG37.5 xb.AG50 xb.AG62.5
## [1,]  1 0.000000 0.0000000000       0         0
## [2,]  1 1.970747 0.0001272736       0         0
## [3,]  1 3.882989 0.0005091021       0         0
## [4,]  1 5.736726 0.0011454777       0         0
## [5,]  1 7.531957 0.0020364083       0         0
## [6,]  1 9.268683 0.0031818860       0         0

U3 <- chol(vcov(n3.qmc))
pred3$se <- sqrt(apply(F3%*%t(U3), 1, function(x) sum(x^2)))
zval <- qnorm(p=c(lwr=0.025, fit=0.5, upr=0.975))
me <- outer(pred3$se, zval, "*")

fit <- predict(n3.qmc, newdata=pred3, level=0)
pred3 <- cbind(pred3, sweep(me, 1, fit, "+"))
str(pred3)

## 'data.frame':    273 obs. of  6 variables:
##  $ x  : num  0 2 4 6 8 10 12 14 16 18 ...
##  $ AG : Factor w/ 3 levels "37.5","50","62.5": 1 1 1 1 1 1 1 1 1 1 ...
##  $ se : num  0.652 0.625 0.601 0.581 0.563 ...
##  $ lwr: num  13.5 14.1 14.7 15.3 15.9 ...
##  $ fit: num  14.7 15.3 15.9 16.4 17 ...
##  $ upr: num  16 16.5 17.1 17.6 18.1 ...

print(p0)+
    as.layer(xyplot(fit~x|AG, data=pred3, type="l", col=2,
                    prepanel=prepanel.cbH, cty="bands",
                    ly=pred3$lwr, uy=pred3$upr,
                    panel=panel.cbH
                    ))

##-----------------------------------------------------------------------------
## Só os valores preditos.

xyplot(fit~x, groups=AG, data=pred1, type="l",
       prepanel=prepanel.cbH, cty="bands", fill=1, alpha=0.2,
       ly=pred1$lwr, uy=pred1$upr,
       panel.groups=panel.cbH, panel=panel.superpose)

xyplot(fit~x, groups=AG, data=pred3, type="l",
       prepanel=prepanel.cbH, cty="bands", fill=1, alpha=0.2,
       ly=pred3$lwr, uy=pred3$upr,
       panel.groups=panel.cbH, panel=panel.superpose)

##-----------------------------------------------------------------------------
## Compara os modelos pelas bandas.

update(p0, strip=strip.custom(bg="gray90"))+
    as.layer(xyplot(fit~x|AG, data=pred1, type="l", col=2,
                    prepanel=prepanel.cbH, cty="bands",
                    ly=pred1$lwr, uy=pred1$upr, fill="red",
                    panel=panel.cbH
                    ))+
    as.layer(xyplot(fit~x|AG, data=pred3, type="l", col=4,
                    prepanel=prepanel.cbH, cty="bands",
                    ly=pred3$lwr, uy=pred3$upr, fill="blue",
                    panel=panel.cbH
                    ))+
    layer(panel.abline(v=c(c1[3], c3[3]), col=c(2,4)), packets=1)+
    layer(panel.abline(v=c(sum(c1[3:4]), c3[4]), col=c(2,4)), packets=2)+
    layer(panel.abline(v=c(sum(c1[c(3,5)]), c3[5]), col=c(2,4)), packets=3)

##-----------------------------------------------------------------------------
## Qual é o melhor modelo afinal?
## Compara as log-verossimilhanças dos modelos finais.

logLik(n1.lpr)

## 'log Lik.' -166.172 (df=7)

logLik(n3.qmr)

## 'log Lik.' -167.8278 (df=7)

AIC(n1.lpr)

##         
## 346.344

AIC(n3.qmr)

##          
## 349.6555

##-----------------------------------------------------------------------------
## Diagnóstico.

r <- residuals(n3.qmr)
f <- fitted(n3.qmr)
plot(r~f)

qqnorm(r)

Parcela subdividida na profundidade

##-----------------------------------------------------------------------------
## Leitura dos dados.

da <- read.table("http://www.leg.ufpr.br/~walmes/data/sistema_gesso_solo.txt",
                 header=TRUE, sep="\t")
## da <- subset(da, sistema=="PD")
da <- subset(da, sistema=="PC")
da <- subset(da, select=-sistema)
str(da)

## 'data.frame':    40 obs. of  17 variables:
##  $ gesso : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ prof  : int  5 5 5 5 10 10 10 10 15 15 ...
##  $ bloco : int  1 2 3 4 1 2 3 4 1 2 ...
##  $ Al    : num  0 0 0 0 0 0 0 1.2 0.1 0 ...
##  $ MO    : num  35.8 30.5 35.5 37.6 35.2 28.2 34.8 36.3 30.5 32.8 ...
##  $ pHKCl : num  5.4 5.2 5.4 4.8 5.5 5.3 5.4 4.5 5.5 5.3 ...
##  $ pHCaCl: num  5.3 5.6 5.6 4.8 5.7 5.8 5.7 4.7 5.8 5.8 ...
##  $ pHH2O : num  6 6.1 6.3 5.6 6.4 6.2 6.3 5.3 6.4 6.3 ...
##  $ Ca    : num  63.9 60.3 65.6 51.3 69.9 62.9 66.6 51.2 58.3 57.6 ...
##  $ Mg    : num  27 25.6 29.3 27.8 32.3 30.3 27.3 25 29.3 30.3 ...
##  $ Hal   : int  34 36 31 55 31 34 33 32 29 31 ...
##  $ S     : num  8.1 3.4 2.6 2.4 8.8 3.3 0.2 1.1 1.3 0.2 ...
##  $ P     : int  34 26 20 18 19 32 22 15 4 10 ...
##  $ K     : num  9.5 10.7 8.7 11.8 3.6 4.2 4.4 4.6 1.8 1.9 ...
##  $ SB    : num  100.4 96.6 103.6 90.9 105.8 ...
##  $ CTC   : num  134 133 135 146 137 ...
##  $ SAT   : int  75 73 77 62 77 74 75 72 76 74 ...

##-----------------------------------------------------------------------------
## Editar.

da <- transform(da,
                bloco=as.factor(bloco),
                gesso=as.factor(gesso),
                Prof=as.factor(prof),
                parc=interaction(bloco, gesso))
str(da)

## 'data.frame':    40 obs. of  19 variables:
##  $ gesso : Factor w/ 2 levels "0","2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ prof  : int  5 5 5 5 10 10 10 10 15 15 ...
##  $ bloco : Factor w/ 4 levels "1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
##  $ Al    : num  0 0 0 0 0 0 0 1.2 0.1 0 ...
##  $ MO    : num  35.8 30.5 35.5 37.6 35.2 28.2 34.8 36.3 30.5 32.8 ...
##  $ pHKCl : num  5.4 5.2 5.4 4.8 5.5 5.3 5.4 4.5 5.5 5.3 ...
##  $ pHCaCl: num  5.3 5.6 5.6 4.8 5.7 5.8 5.7 4.7 5.8 5.8 ...
##  $ pHH2O : num  6 6.1 6.3 5.6 6.4 6.2 6.3 5.3 6.4 6.3 ...
##  $ Ca    : num  63.9 60.3 65.6 51.3 69.9 62.9 66.6 51.2 58.3 57.6 ...
##  $ Mg    : num  27 25.6 29.3 27.8 32.3 30.3 27.3 25 29.3 30.3 ...
##  $ Hal   : int  34 36 31 55 31 34 33 32 29 31 ...
##  $ S     : num  8.1 3.4 2.6 2.4 8.8 3.3 0.2 1.1 1.3 0.2 ...
##  $ P     : int  34 26 20 18 19 32 22 15 4 10 ...
##  $ K     : num  9.5 10.7 8.7 11.8 3.6 4.2 4.4 4.6 1.8 1.9 ...
##  $ SB    : num  100.4 96.6 103.6 90.9 105.8 ...
##  $ CTC   : num  134 133 135 146 137 ...
##  $ SAT   : int  75 73 77 62 77 74 75 72 76 74 ...
##  $ Prof  : Factor w/ 5 levels "5","10","15",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ parc  : Factor w/ 8 levels "1.0","2.0","3.0",..: 1 2 3 4 1 2 3 4 1 2 ...

##-----------------------------------------------------------------------------
## Ver.

da$y <- da$P

xyplot(y~prof, groups=gesso, data=da, type=c("p","a"),
       auto.key=TRUE, jitter.x=TRUE)

xyplot(y~prof|gesso, groups=parc, data=da, type=c("p","a"),
       auto.key=TRUE, jitter.x=TRUE)

##-----------------------------------------------------------------------------
## Ajuste do modelo de parcela subdivididas (assume independência entre
## observações de uma mesma parcela na profundidade).

m0 <- lme(y~bloco+gesso*Prof,
          random=~1|parc, data=da,
          method="ML")

##-----------------------------------------------------------------------------
## Diagnóstico.

r <- residuals(m0, type="pearson")
f <- fitted(m0)
bp <- unlist(ranef(m0)[1])

grid.arrange(ncol=2,
             xyplot(r~f, type=c("p", "smooth")),
             xyplot(sqrt(abs(r))~f, type=c("p", "smooth")),
             qqmath(r), qqmath(bp))

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

## Embora esteja considerando só a parte de efeito aleatório um nível
## anterior ao de resíduo, serve de indicação.
MASS::boxcox(lm(y~parc, data=da))

##-----------------------------------------------------------------------------
## Modelo com a variável transformada. Refazer a análise dos resíduos.

m0 <- update(m0, log(y)~.)

r <- residuals(m0, type="pearson")
f <- fitted(m0)
bp <- unlist(ranef(m0)[1])

grid.arrange(ncol=2,
             xyplot(r~f, type=c("p", "smooth")),
             xyplot(sqrt(abs(r))~f, type=c("p", "smooth")),
             qqmath(r), qqmath(bp))

## Teste de Wald para os termos de efeito fixo.
anova(m0)

##             numDF denDF   F-value p-value
## (Intercept)     1    24 1044.9454  <.0001
## bloco           3     3    7.0616  0.0713
## gesso           1     3    0.6589  0.4764
## Prof            4    24  114.8043  <.0001
## gesso:Prof      4    24    3.5606  0.0204

##-----------------------------------------------------------------------------
## Mas e a correlação serial ao longo do perfil?

plot(ACF(m0))

S <- split(data.frame(r=r, prof=da$prof), da$parc)
S <- lapply(S,
            function(i){
                i$r[order(i$prof)]
            })
S <- do.call(rbind, S)

splom(S, type=c("p","r"))

round(cor(S), 2)

##       [,1]  [,2]  [,3]  [,4]  [,5]
## [1,]  1.00 -0.11 -0.69 -0.87  0.52
## [2,] -0.11  1.00  0.59 -0.27 -0.17
## [3,] -0.69  0.59  1.00  0.48 -0.31
## [4,] -0.87 -0.27  0.48  1.00 -0.26
## [5,]  0.52 -0.17 -0.31 -0.26  1.00

## ## Estrutura de correlação não estruturada (tem k*(k-1)/2 parâmetros).
## m1 <- update(m0,
##              correlation=corSymm(form=~as.integer(Prof)|parc))
## summary(m1)

## anova(m0, m1)

## Estrutura de continous AR(1).
m2 <- update(m0, correlation=corCAR1(form=~prof|parc))
anova(m0, m2)

##    Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## m0     1 15 38.01346 63.34665 -4.006729                        
## m2     2 16 38.01957 65.04164 -3.009786 1 vs 2 1.993886  0.1579

## AIC(m1)
AIC(m2)

## [1] 38.01957

##-----------------------------------------------------------------------------
## Os resultados podem mudar.

anova(m0)

##             numDF denDF   F-value p-value
## (Intercept)     1    24 1044.9454  <.0001
## bloco           3     3    7.0616  0.0713
## gesso           1     3    0.6589  0.4764
## Prof            4    24  114.8043  <.0001
## gesso:Prof      4    24    3.5606  0.0204

## anova(m1)
anova(m2)

##             numDF denDF  F-value p-value
## (Intercept)     1    24 759.6583  <.0001
## bloco           3     3   5.2934  0.1022
## gesso           1     3   0.4484  0.5510
## Prof            4    24  98.5976  <.0001
## gesso:Prof      4    24   3.8457  0.0149

##-----------------------------------------------------------------------------
## Qual foi a correlação estimada?

## summary(m1)$modelStruct$corStruct
summary(m2)$modelStruct$corStruct

## Correlation structure of class corCAR1 representing
##       Phi 
## 0.7373746

##-----------------------------------------------------------------------------
## Obter valores preditos para profundidade em cada nível de gesso e
## sistema.

X <- LSmatrix(lm(formula(m2), data=da),
              effect=c("gesso","Prof"))

grid <- attr(X, "grid")

## Médias com IC.
ci <- confint(glht(m2, linfct=X), calpha=univariate_calpha())

grid <- cbind(grid, ci$confint)
grid <- equallevels(grid, da)
grid

##    gesso Prof  Estimate       lwr       upr
## 1      0    5 3.1676403 2.9069547 3.4283259
## 2      2    5 3.0485000 2.7878144 3.3091856
## 3      0   10 3.0523169 2.7916313 3.3130025
## 4      2   10 2.4393263 2.1786407 2.7000118
## 5      0   15 1.6153670 1.3546815 1.8760526
## 6      2   15 1.3540251 1.0933395 1.6147107
## 7      0   20 0.3465736 0.0858880 0.6072592
## 8      2   20 0.9222199 0.6615343 1.1829055
## 9      0   25 0.3465736 0.0858880 0.6072592
## 10     2   25 0.3465736 0.0858880 0.6072592

##-----------------------------------------------------------------------------
## Representação.

## l <- sapply(levels(grid$gesso),
##             function(l){
##                 eval(substitute(expression(gesso~(ton~ha^{-1})),
##                                 list(gesso=l)))
##             })

l <- levels(grid$gesso)
key <- list(type="o", divide=1, columns=length(l),
            title=expression("Gesso agrícola"~(ton~ha^{-1})),
            cex.title=1.1,
            lines=list(pch=1:length(l),
                col=trellis.par.get("plot.symbol")$col),
            text=list(l))

## Para que os pontos sejam corretamente representados, deve-se
## reordenar os dados.
grid <- grid[with(grid, order(Prof, gesso)),]

segplot(Prof~lwr+upr, horizontal=FALSE,
        centers=Estimate, groups=grid$gesso, f=0.05, layout=c(NA,1),
        data=grid, draw=FALSE, panel=panel.segplotBy, key=key,
        pch=as.integer(grid$gesso),
        ylab="log Fósforo",
        xlab="Camada do solo (cm)")

O R só falta falar, não é?

##-----------------------------------------------------------------------------
## Fazendo o R falar.

source("~/Dropbox/CursoR/cpaa/sayHello.R")
sayHello()
sayHello("Estou adorando estar aqui com vocês.")

st <- Sys.time()
fala <- sprintf("Agora são %s horas e %s minutos. A sua análise acaba de ser concluída",
        format(st, "%H"), format(st, "%M"))
sayHello(fala)

Veja um exemplo de análise de modelo polinomial com termos de efeito aleatório aqui: http://www.leg.ufpr.br/~walmes/analises/RHTBGoes/copaiba.html .

##-----------------------------------------------------------------------------
## Informações da sessão.

Sys.time()

## [1] "2014-12-11 21:05:15 BRST"

sessionInfo()

## R version 3.1.2 (2014-10-31)
## 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   grid      stats     graphics  grDevices utils     datasets 
## [9] base     
## 
## other attached packages:
##  [1] wzRfun_0.4          rootSolve_1.6.5.1   nlme_3.1-118        multcomp_1.3-7     
##  [5] TH.data_1.0-5       mvtnorm_1.0-1       doBy_4.5-12         survival_2.37-7    
##  [9] reshape_0.8.5       plyr_1.8.1          alr3_2.0.5          car_2.0-22         
## [13] gridExtra_0.9.1     latticeExtra_0.6-26 RColorBrewer_1.0-5  lattice_0.20-29    
## [17] rmarkdown_0.3.3     knitr_1.8          
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.4    evaluate_0.5.5  formatR_1.0     htmltools_0.2.6 MASS_7.3-35    
##  [6] Matrix_1.1-4    nnet_7.3-8      Rcpp_0.11.3     sandwich_2.3-2  stringr_0.6.2  
## [11] tools_3.1.2     yaml_2.1.13     zoo_1.7-11

Aplicação de modelos de regressão linear e não linear em ciências agrárias
09 à 11 de Dezembro de 2014 - Goiânia - GO	Prof. Dr. Walmes M. Zeviani
Embrapa Arroz e Feijão	Lab. de Estatística e Geoinformação - LEG
	Departamento de Estatística - UFPR

Aplicação de modelos de regressão linear e não linear em ciências agrárias

Regressão na presença de termos de efeito aleatório

Fatorial duplo analisado com regressão não linear (misto)

Parcela subdividida na profundidade

O R só falta falar, não é?