1019inclass1

1 load package

pacman::p_load(tidyverse, nlme)

2 data input and management

# data input
dta <- read.csv("ade.csv")

# see the data structure
str(dta)

## 'data.frame':    60 obs. of  3 variables:
##  $ Dog : chr  "D11" "D11" "D11" "D11" ...
##  $ Hour: num  0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 ...
##  $ ADE : num  5.7 4.7 4.8 4.9 3.88 3.51 2.8 2.6 2.5 2.5 ...

# show the first six data
head(dta)

##   Dog Hour  ADE
## 1 D11  0.0 5.70
## 2 D11  0.5 4.70
## 3 D11  1.0 4.80
## 4 D11  1.5 4.90
## 5 D11  2.0 3.88
## 6 D11  2.5 3.51

#產生hour的factor變項
dta <- mutate(dta, Hour_f = factor(Hour))

# 設定佈景主題
theme_set(theme_bw())

3 means with error bars

ggplot(dta, aes(Hour, ADE)) +
 geom_point(alpha = .3) +
 stat_summary(fun.y = mean, geom = "line") +
 stat_summary(fun.y = mean, geom = "point") +
 stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
 labs(x = "Time (in hours)", y = "ADE")

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: `fun.y` is deprecated. Use `fun` instead.

4 model

#constant variance function structure
# correlation=Defaults to NULL, corresponding to uncorrelated errors.
summary(m1 <- gls(ADE ~ Hour,
                  weights = varIdent(form = ~ 1 | Hour_f),
                  data = dta))

## Generalized least squares fit by REML
##   Model: ADE ~ Hour 
##   Data: dta 
##        AIC      BIC    logLik
##   196.6523 221.3776 -86.32617
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | Hour_f 
##  Parameter estimates:
##          0        0.5          1        1.5          2        2.5          3 
## 1.00000000 0.55447427 0.64132902 0.51921510 0.43187908 0.08245758 0.10901472 
##        3.5          4        4.5 
## 0.17046144 0.21770576 0.11133075 
## 
## Coefficients:
##                 Value  Std.Error   t-value p-value
## (Intercept)  5.435552 0.27441443 19.807821       0
## Hour        -0.701874 0.08386437 -8.369153       0
## 
##  Correlation: 
##      (Intr)
## Hour -0.964
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.3854719 -0.6272929 -0.1826742  0.5824007  2.3554705 
## 
## Residual standard error: 3.508619 
## Degrees of freedom: 60 total; 58 residual

# corStruct = corAR1
summary(m2 <- gls(ADE ~ Hour,
                   weights = varIdent(form = ~ 1 | Hour_f),
                   correlation = corAR1(form = ~ Hour | Dog),
                   data = dta))

## Generalized least squares fit by REML
##   Model: ADE ~ Hour 
##   Data: dta 
##        AIC      BIC    logLik
##   168.5913 195.3771 -71.29566
## 
## Correlation Structure: ARMA(1,0)
##  Formula: ~Hour | Dog 
##  Parameter estimate(s):
##      Phi1 
## 0.7553498 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | Hour_f 
##  Parameter estimates:
##          0        0.5          1        1.5          2        2.5          3 
## 1.00000000 0.32803241 0.37194739 0.28590639 0.34018814 0.12571538 0.09896982 
##        3.5          4        4.5 
## 0.16894933 0.22219099 0.09978382 
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept)  4.749848 0.3999065 11.877395       0
## Hour        -0.599160 0.1096671 -5.463442       0
## 
##  Correlation: 
##      (Intr)
## Hour -0.94 
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.43566691 -0.09530459  0.29216375  1.34205811  5.06634565 
## 
## Residual standard error: 3.278488 
## Degrees of freedom: 60 total; 58 residual

# corStruct = corCompSymm
summary(m3 <- gls(ADE ~ Hour,
                  weights = varIdent(form = ~ 1 | Hour_f),
                  correlation = corCompSymm(form = ~ 1 | Dog),
                  data = dta))

## Generalized least squares fit by REML
##   Model: ADE ~ Hour 
##   Data: dta 
##        AIC      BIC    logLik
##   161.6077 188.3935 -67.80387
## 
## Correlation Structure: Compound symmetry
##  Formula: ~1 | Dog 
##  Parameter estimate(s):
##       Rho 
## 0.7217933 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | Hour_f 
##  Parameter estimates:
##         0       0.5         1       1.5         2       2.5         3       3.5 
## 1.0000000 0.5100528 0.6140381 0.5228105 0.4479331 0.1852681 0.1129502 0.1611728 
##         4       4.5 
## 0.2309740 0.1131425 
## 
## Coefficients:
##                 Value  Std.Error   t-value p-value
## (Intercept)  4.494975 0.30766482 14.609974       0
## Hour        -0.529169 0.06981253 -7.579863       0
## 
##  Correlation: 
##      (Intr)
## Hour -0.957
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -0.66528429 -0.02390899  0.26405734  0.89073611  2.44460225 
## 
## Residual standard error: 3.832969 
## Degrees of freedom: 60 total; 58 residual

# exponential variance function structure
summary(m3a <- gls(ADE ~ Hour,
                  weights = varExp(form = ~ Hour),
                  correlation = corCompSymm(form = ~ 1 | Dog),
                  data = dta))

## Generalized least squares fit by REML
##   Model: ADE ~ Hour 
##   Data: dta 
##        AIC      BIC    logLik
##   163.9703 174.2725 -76.98516
## 
## Correlation Structure: Compound symmetry
##  Formula: ~1 | Dog 
##  Parameter estimate(s):
##       Rho 
## 0.6145863 
## Variance function:
##  Structure: Exponential of variance covariate
##  Formula: ~Hour 
##  Parameter estimates:
##      expon 
## -0.4204529 
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept)  5.302798 0.6777356  7.824287       0
## Hour        -0.674162 0.1346238 -5.007751       0
## 
##  Correlation: 
##      (Intr)
## Hour -0.982
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.1017685 -0.3945215 -0.1031040  0.3711706  2.7263751 
## 
## Residual standard error: 3.079988 
## Degrees of freedom: 60 total; 58 residual

# random component and constant variance function structure           
summary(m3b <- lme(ADE ~ Hour, random = ~ 1 | Dog,
                  weights = varIdent(form = ~ 1 | Hour_f),
                  correlation = corCompSymm(form = ~ 1 | Dog),
                  data = dta))

## Linear mixed-effects model fit by REML
##  Data: dta 
##        AIC     BIC    logLik
##   153.2468 182.093 -62.62338
## 
## Random effects:
##  Formula: ~1 | Dog
##         (Intercept) Residual
## StdDev:   0.3654553 3.928157
## 
## Correlation Structure: Compound symmetry
##  Formula: ~1 | Dog 
##  Parameter estimate(s):
##       Rho 
## 0.8457657 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | Hour_f 
##  Parameter estimates:
##          0        0.5          1        1.5          2        2.5          3 
## 1.00000000 0.54618337 0.63785632 0.53997928 0.47199948 0.16064282 0.06300446 
##        3.5          4        4.5 
## 0.21603493 0.25337758 0.18268767 
## Fixed effects: ADE ~ Hour 
##                 Value  Std.Error DF   t-value p-value
## (Intercept)  4.811462 0.23464110 53  20.50562       0
## Hour        -0.612618 0.05378694 53 -11.38973       0
##  Correlation: 
##      (Intr)
## Hour -0.745
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -0.9223391 -0.1299528  0.1645922  0.6995419  2.4095565 
## 
## Number of Observations: 60
## Number of Groups: 6

5 comparing models

#m1=constant variance function and unstructured,m2=ARstructured
anova(m1, m2)

##    Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## m1     1 12 196.6523 221.3776 -86.32617                        
## m2     2 13 168.5913 195.3771 -71.29566 1 vs 2 30.06102  <.0001

#m3=constant variance function and compound structured, m3a=exponential variance function and compound structured
anova(m1, m3, m3a)

##     Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## m1      1 12 196.6523 221.3776 -86.32617                        
## m3      2 13 161.6077 188.3935 -67.80387 1 vs 2 37.04460  <.0001
## m3a     3  5 163.9703 174.2725 -76.98516 2 vs 3 18.36258  0.0187

6 residual plot

plot(m3, resid(., type = "pearson") ~ fitted(.) | Dog,
     layout = c(6, 1), aspect = 2, abline = 0, pch = 20, cex = .8,
     xlab = "Ftted values", ylab = "Pearson residuals")

7 The end