Week7 In-class exercises 2
Ching-Fang Wu
Sun Oct 25 17:22:59 2020
The arrhythmogenic dose of epinephrine (ADE) is determined by infusing epinephrine into an animal until arrhythmia criterion is reached. The arrhythmia criterion was the occurrence of four intermittent or continuous premature ventricular contractions. Once the criterion has been reached, the infusion required to produce it is recorded. The ADE is then calculated by multiplying the duration of the infusion by the infusion rate. The data give the ADE for six dogs measured at ten consecutive times.
Source: Shoukri, M.M., & Chaudhary, M.A. (2007). Analysis of Correlated Data with SAS and R. p. 218.
Column 1: Dog ID Column 2: Time in hours Column 3: Arrhythmogenic dose of epinephrine (ADE)
1 Data management
pacman::p_load(tidyverse, nlme) #pacman::p_load()可同時取代 install.packages("A","B") 和 library(A, B)#input data
dta <- read.csv("C:/Users/Ching-Fang Wu/Documents/data/ade.csv")#show strurcture of data
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 first and tail 6 lines
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.51tail(dta)## Dog Hour ADE
## 55 D16 2.0 3.9
## 56 D16 2.5 4.0
## 57 D16 3.0 3.0
## 58 D16 3.5 2.4
## 59 D16 4.0 2.1
## 60 D16 4.5 2.0dta <- mutate(dta, Hour_f = factor(Hour)) #新增時間類別變數Hour_ftheme_set(theme_bw()) #2 Data Visualization
2.1 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.
mean及standard deviation隨著時間遞減。
3 Modeling
3.1 m1
#m1考量每小時的不同變異
summary(m1 <- gls(ADE ~ Hour,
weights = varIdent(form = ~ 1 | Hour_f),
data = dta)) #varIdent代表Different variances per stratum，form = ~ 1 | Hour_f 的意思是Constant variance within stratum.## 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 residual3.2 m2
summary(m2 <- gls(ADE ~ Hour,
weights = varIdent(form = ~ 1 | Hour_f),
correlation = corAR1(form = ~ Hour | Dog), #這裡的共變異數矩陣為一階自我相關corAR1，每筆資料和上一個量測有關。此外，以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 residual3.3 m3
summary(m3 <- gls(ADE ~ Hour,
weights = varIdent(form = ~ 1 | Hour_f),
correlation = corCompSymm(form = ~ 1 | Dog),
data = dta)) #這裡的共變異數矩陣為Compound Symmetry Correlation Structure，假設重複測量之間的相關性都一樣。## 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 residual3.4 m3a
summary(m3a <- gls(ADE ~ Hour,
weights = varExp(form = ~ Hour), #varExp是Exponential variance structure，varExp(form = ~ Hour)表示變異數隨時間而不同
correlation = corCompSymm(form = ~ 1 | Dog), #corCompSymm假設重複測量之間的相關性都一樣
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 residual3.5 m3b
#mixed effect model
summary(m3b <- lme(ADE ~ Hour, random = ~ 1 | Dog, #random effect
weights = varIdent(form = ~ 1 | Hour_f),
correlation = corCompSymm(form = ~ 1 | Dog),##corCompSymm假設重複測量之間的相關性都一樣
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##comparing models
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比較m1和m2，選AIC值最小值m2。
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比較m1、m3和m3a，選AIC最小值：m3
4 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")
#The end