Carregando Pacotes

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

library(lattice)
library(latticeExtra)
library(splines) # bs(), ns().
library(mgcv)    # gam()
library(wzRfun)  # panel.cbH().

Carregando os Dados

Os são referentes ao monitoramento da temperatura do corpo de suínos em um ambiente preabate com aspersão de água e sem aspersão. A aspersão visa dimunir a temperatura do animal e era feita por 30 minutos a cada 30 minutos.

#-----------------------------------------------------------------------
# Leitura.

sui <- read.table("http://www.leg.ufpr.br/~walmes/data/preabate.txt",
                  header = TRUE, dec = ",")
sui$trat <- factor(-1 * sui$trat)
levels(sui$trat) <- c("Sem aspersão", "Com aspersão")
str(sui)
## 'data.frame':    384 obs. of  5 variables:
##  $ trat: Factor w/ 2 levels "Sem aspersão",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ asp : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ hora: int  0 0 0 0 5 5 5 5 10 10 ...
##  $ rep : int  1 2 3 4 1 2 3 4 1 2 ...
##  $ temp: num  29 30.9 34.5 31.3 29.6 ...
xtabs(~trat + rep, data = sui)
##               rep
## trat            1  2  3  4
##   Sem aspersão 48 48 48 48
##   Com aspersão 48 48 48 48
xyplot(temp ~ hora,
       groups = trat,
       data = sui,
       xlab = "Minutos à partir das 07:00 da manhã",
       ylab = "Temperatura corporal (ºC)",
       auto.key = TRUE) +
    glayer(panel.smoother(..., span = 0.4))

xyplot(temp ~ hora | trat, groups = asp, data = sui,
       type = c("p", "smooth"))

# Efeito de animal?
xyplot(temp ~ hora | rep + trat,
       data = sui,
       type = "b")

xyplot(temp ~ hora | trat,
       groups = rep,
       data = sui,
       type = "b")

Regressão Local

#-----------------------------------------------------------------------
# Loess.

# Apenas com os dados do tratamento com aspersão.
ca <- subset(sui, trat == "Com aspersão")

# m0 <- loess(temp ~ hora, data = ca)
# m0 <- loess(temp ~ hora, data = ca, span = 0.25, degree = 1)
m0 <- loess(temp ~ hora, data = ca, span = 0.25)

# Medidas de ajuste.
summary(m0)
## Call:
## loess(formula = temp ~ hora, data = ca, span = 0.25)
## 
## Number of Observations: 192 
## Equivalent Number of Parameters: 11.87 
## Residual Standard Error: 1.307 
## Trace of smoother matrix: 13.12  (exact)
## 
## Control settings:
##   span     :  0.25 
##   degree   :  2 
##   family   :  gaussian
##   surface  :  interpolate      cell = 0.2
##   normalize:  TRUE
##  parametric:  FALSE
## drop.square:  FALSE
# Classe e métodos para a classe.
class(m0)
## [1] "loess"
methods(class = class(m0))
## [1] anova   predict print   summary
## see '?methods' for accessing help and source code
qqnorm(residuals(m0))

# Predição.
pred <- data.frame(hora = seq(min(ca$hora), max(ca$hora), l = 200))
aux <- predict(m0, newdata = pred, se = TRUE)
aux$me <- outer(aux$se.fit,
                qnorm(0.975) * c(lwr = -1, fit = 0, upr = 1),
                FUN = "*")

aux <- sweep(aux$me, MARGIN = 1, STATS = aux$fit, FUN = "+")
pred <- cbind(pred, aux)
str(pred)
## 'data.frame':    200 obs. of  4 variables:
##  $ hora: num  0 1.18 2.36 3.54 4.72 ...
##  $ lwr : num  30.5 30.2 30 29.8 29.5 ...
##  $ fit : num  31.5 31.2 30.8 30.5 30.2 ...
##  $ upr : num  32.5 32.1 31.7 31.3 30.9 ...
p1 <- xyplot(temp ~ hora,
             data = ca,
             ylab = "Temperatura do corpo",
             xlab = "Minutos a partir das 7h00")

p2 <- xyplot(fit ~ hora,
             data = pred,
             type = "l",
             ly = pred$lwr,
             uy = pred$upr,
             cty = "bands",
             prepanel = prepanel.cbH,
             panel = panel.cbH)

p1 + as.layer(p2)

Splines

#-----------------------------------------------------------------------
# Splines.

# m0 <- lm(temp ~ bs(hora), data = ca)
# m0 <- lm(temp ~ bs(hora, degree = 13, df = 1), data = ca)
m0 <- lm(temp ~ ns(hora, df = 13), data = ca)
summary(m0)
## 
## Call:
## lm(formula = temp ~ ns(hora, df = 13), data = ca)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7022 -0.9033 -0.0045  0.7134  3.8193 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         31.66367    0.55157  57.406  < 2e-16 ***
## ns(hora, df = 13)1  -2.27264    0.85325  -2.664 0.008442 ** 
## ns(hora, df = 13)2  -1.22015    0.97361  -1.253 0.211771    
## ns(hora, df = 13)3  -3.37650    0.89413  -3.776 0.000217 ***
## ns(hora, df = 13)4  -1.89242    0.98766  -1.916 0.056959 .  
## ns(hora, df = 13)5   1.24033    0.91510   1.355 0.177005    
## ns(hora, df = 13)6  -2.62154    0.92759  -2.826 0.005249 ** 
## ns(hora, df = 13)7  -3.10635    0.96947  -3.204 0.001605 ** 
## ns(hora, df = 13)8   2.75124    0.92410   2.977 0.003314 ** 
## ns(hora, df = 13)9   0.08984    0.92098   0.098 0.922400    
## ns(hora, df = 13)10 -4.65297    0.96546  -4.819 3.07e-06 ***
## ns(hora, df = 13)11  1.10806    0.80262   1.381 0.169149    
## ns(hora, df = 13)12 -3.62265    1.42213  -2.547 0.011701 *  
## ns(hora, df = 13)13  3.78342    0.66592   5.682 5.36e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.304 on 178 degrees of freedom
## Multiple R-squared:  0.5533, Adjusted R-squared:  0.5207 
## F-statistic: 16.96 on 13 and 178 DF,  p-value: < 2.2e-16
# Onde foram posicionados os nós.
n <- c(attr(m0$model[[2]], "knots"),
       attr(m0$model[[2]], "Boundary.knots"))
cbind(n)
##             n
## 7.692308%  15
## 15.38462%  35
## 23.07692%  55
## 30.76923%  70
## 38.46154%  90
## 46.15385% 110
## 53.84615% 125
## 61.53846% 145
## 69.23077% 165
## 76.92308% 180
## 84.61538% 200
## 92.30769% 220
##             0
##           235
# Predição.
pred <- data.frame(hora = seq(min(ca$hora), max(ca$hora), l = 200))
aux <- predict(m0, newdata = pred, interval = "confidence")
pred <- cbind(pred, aux)
str(pred)
## 'data.frame':    200 obs. of  4 variables:
##  $ hora: num  0 1.18 2.36 3.54 4.72 ...
##  $ fit : num  31.7 31.3 31 30.6 30.3 ...
##  $ lwr : num  30.6 30.3 30.1 29.8 29.6 ...
##  $ upr : num  32.8 32.3 31.8 31.4 31 ...
p1 <- xyplot(temp ~ hora,
             data = ca,
             ylab = "Temperatura do corpo",
             xlab = "Minutos a partir das 7h00")
p2 <- xyplot(fit ~ hora,
             data = pred,
             type = "l",
             ly = pred$lwr,
             uy = pred$upr,
             cty = "bands",
             prepanel = prepanel.cbH,
             panel = panel.cbH)
p1 + as.layer(p2) +
    layer(panel.abline(v = n, lty = 2, col = "gray50"))

# Usando um knots a cada 30 minutos.
knots <- seq(min(ca$hora), max(ca$hora), by = 30)[-1]
length(knots)
## [1] 7
m0 <- lm(temp ~ ns(hora, knots = knots), data = ca)
summary(m0)
## 
## Call:
## lm(formula = temp ~ ns(hora, knots = knots), data = ca)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2204 -1.1042 -0.2266  1.0065  4.3262 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               30.3920     0.5490  55.361  < 2e-16 ***
## ns(hora, knots = knots)1  -0.6880     0.7682  -0.896    0.372    
## ns(hora, knots = knots)2  -0.1575     0.9505  -0.166    0.869    
## ns(hora, knots = knots)3   0.3780     0.8646   0.437    0.662    
## ns(hora, knots = knots)4  -0.6855     0.9164  -0.748    0.455    
## ns(hora, knots = knots)5   3.7409     0.8993   4.160 4.89e-05 ***
## ns(hora, knots = knots)6  -3.2224     0.7710  -4.179 4.52e-05 ***
## ns(hora, knots = knots)7  -1.0871     1.4165  -0.767    0.444    
## ns(hora, knots = knots)8   4.7919     0.6642   7.215 1.39e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.558 on 183 degrees of freedom
## Multiple R-squared:  0.3447, Adjusted R-squared:  0.3161 
## F-statistic: 12.03 on 8 and 183 DF,  p-value: 9.489e-14

Modelos aditivos generalizados (GAM)

#-----------------------------------------------------------------------
# Ajustar GAM.

m0 <- gam(temp ~ s(hora), data = ca)
summary(m0)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## temp ~ s(hora)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  30.2432     0.0951     318   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df     F p-value    
## s(hora) 8.92  8.998 22.55  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.511   Deviance explained = 53.3%
## GCV = 1.8311  Scale est. = 1.7365    n = 192
pred <- data.frame(hora = seq(min(ca$hora), max(ca$hora), l = 200))
aux <- predict(m0, newdata = pred, se = TRUE)
aux$me <- outer(aux$se.fit,
                c(lwr = -1, fit = 0, upr = 1) * qnorm(0.975),
                FUN = "*")
aux <- sweep(aux$me, MARGIN = 1, STATS = aux$fit, FUN = "+")
pred <- cbind(pred, aux)
str(pred)
## 'data.frame':    200 obs. of  4 variables:
##  $ hora: num  0 1.18 2.36 3.54 4.72 ...
##  $ lwr : num  30.5 30.3 30.1 29.8 29.6 ...
##  $ fit : num  31.5 31.2 30.9 30.6 30.3 ...
##  $ upr : num  32.5 32.1 31.7 31.4 31 ...
p1 <- xyplot(temp ~ hora,
             data = ca,
             ylab = "Temperatura do corpo",
             xlab = "Minutos a partir das 7h00")
p2 <- xyplot(fit ~ hora,
             data = pred,
             type = "l",
             ly = pred$lwr,
             uy = pred$upr,
             cty = "bands",
             prepanel = prepanel.cbH,
             panel = panel.cbH)
p1 + as.layer(p2)

#-----------------------------------------------------------------------
# Ajuste para os dois níveis de tratamento.

# Modelo sem interação trat com hora.
g1 <- gam(temp ~ trat + s(hora), data = sui)

# Modelo em que em cada trat a hora tem um efeito diferente.
g0 <- gam(temp ~ trat + s(hora, by = trat), data = sui)

# Para hipótese de interação nula.
anova(g0, g1, test = "F")
## Analysis of Deviance Table
## 
## Model 1: temp ~ trat + s(hora, by = trat)
## Model 2: temp ~ trat + s(hora)
##   Resid. Df Resid. Dev      Df Deviance      F    Pr(>F)    
## 1    367.71     562.50                                      
## 2    373.01     824.09 -5.3052  -261.59 32.326 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Resumo do ajuste.
summary(g0)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## temp ~ trat + s(hora, by = trat)
## 
## Parametric coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      34.76667    0.08913  390.06   <2e-16 ***
## tratCom aspersão -4.52344    0.12605  -35.89   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                            edf Ref.df     F p-value    
## s(hora):tratSem aspersão 4.302  5.295 92.57  <2e-16 ***
## s(hora):tratCom aspersão 8.930  8.999 25.76  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.84   Deviance explained = 84.5%
## GCV = 1.5883  Scale est. = 1.5253    n = 384
pred <- expand.grid(hora = seq(min(sui$hora), max(sui$hora), l = 100),
                    trat = levels(sui$trat))
aux <- predict(g0, newdata = pred, se = TRUE)
aux$me <- outer(aux$se.fit,
                c(lwr = -1, fit = 0, upr = 1) * qnorm(0.975),
                FUN = "*")
aux <- sweep(aux$me, MARGIN = 1, STATS = aux$fit, FUN = "+")
pred <- cbind(pred, aux)
str(pred)
## 'data.frame':    200 obs. of  5 variables:
##  $ hora: num  0 2.37 4.75 7.12 9.49 ...
##  $ trat: Factor w/ 2 levels "Sem aspersão",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ lwr : num  30.3 30.4 30.6 30.7 30.8 ...
##  $ fit : num  31 31 31.1 31.2 31.3 ...
##  $ upr : num  31.6 31.7 31.7 31.7 31.8 ...
xmn <- seq(0, 4 * 60, by = 30)
xhr <- seq(as.POSIXct("07:00", format = "%H:%M"),
           as.POSIXct("11:00", format = "%H:%M"), "30 min")
format(xhr, "%H:%M")
## [1] "07:00" "07:30" "08:00" "08:30" "09:00" "09:30" "10:00" "10:30"
## [9] "11:00"
p1 <- xyplot(temp ~ hora, groups = trat, data = sui,
             xlab = "Minutos à partir das 07:00 da manhã",
             ylab = "Temperatura corporal (ºC)",
             scales = list(x = list(at = xmn,
                                    labels = format(xhr, "%H:%M"))),
             auto.key = list(columns = 2), panel = function(...) {
                 panel.abline(v = (0:9) * 30, lty = 3, col = "gray50")
                 panel.rect((0:3) * 60, 22, ((0:3) + 0.5) * 60, 40,
                            col = "gray35", border = "transparent",
                            alpha = 0.1)
                 panel.xyplot(...)
             })

p2 <- xyplot(fit ~ hora,
             groups = trat,
             data = pred,
             type = "l",
             prepanel = prepanel.cbH,
             cty = "bands",
             ly = pred$lwr,
             uy = pred$upr,
             alpha = 0.1,
             fill = 1,
             panel.groups = panel.cbH,
             panel = panel.superpose)

p1 + as.layer(p2)

Informações da Sessão

## Updated in 2016-10-10.
## 
## R version 3.3.1 (2016-06-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.1 LTS
## 
## locale:
##  [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=pt_BR.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=pt_BR.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=pt_BR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] splines   stats     graphics  grDevices utils     datasets 
## [7] base     
## 
## other attached packages:
## [1] wzRfun_0.70         mgcv_1.8-13         nlme_3.1-128       
## [4] latticeExtra_0.6-28 RColorBrewer_1.1-2  lattice_0.20-33    
## [7] rmarkdown_1.0       knitr_1.14         
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.7      magrittr_1.5     MASS_7.3-45     
##  [4] multcomp_1.4-6   stringr_1.0.0    tools_3.3.1     
##  [7] rpanel_1.1-3     grid_3.3.1       TH.data_1.0-7   
## [10] htmltools_0.3.5  survival_2.39-4  yaml_2.1.13     
## [13] digest_0.6.10    Matrix_1.2-6     formatR_1.4     
## [16] codetools_0.2-14 evaluate_0.9     sandwich_2.3-4  
## [19] doBy_4.5-15      stringi_1.1.1    methods_3.3.1   
## [22] mvtnorm_1.0-5    zoo_1.7-13