Routine for PH 628
Lecture 13-14: Longitudinal Data Analysis
Gabriel Carrasco Escobar
UC San Diego; Universidad Peruana Cayetano Heredia (UPCH)Inquiries regarding this code contact gabriel.carrasco@upch.pe
rm(list=ls())
library(dplyr); library(ggplot2); library(epiDisplay); library(tidyr); library(haven)
library(GGally); library(afex); library(car); library(forcats); library(nlme)Random coefficients (intercept only) model
Go to SAS output, pp 143-147, 148-153
head(dat)## # A tibble: 6 x 8
## PATIENT UV0 UV4 UV8 GROUP LNUV0 LNUV4 LNUV8
## <chr> <dbl> <dbl> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 M01 2566. 2464. 2708 control 7.85 7.81 7.90
## 2 M02 3468. 3478. NA treatment 8.15 8.15 NA
## 3 M03 2224. 1499 1364. treatment 7.71 7.31 7.22
## 4 M04 1050. 900 1249 treatment 6.96 6.80 7.13
## 5 M05 3616. 4067 3881 control 8.19 8.31 8.26
## 6 M06 2932 2102. 2612. treatment 7.98 7.65 7.87dt<-dat %>% dplyr::select(PATIENT, GROUP, LNUV0, LNUV4, LNUV8) %>% gather(key = TIME, value = LNUV, LNUV0:LNUV8) %>%
mutate_at(c("PATIENT","GROUP","TIME"), funs(as.factor))The MEANS procedure
(des <- dt %>% group_by(GROUP, TIME) %>%
summarize(Obs=n(), N=sum(!is.na(LNUV)), Mean=mean(LNUV,na.rm = T), Std_Dev=sd(LNUV,na.rm = T), Minimum=min(LNUV,na.rm = T), Maximum=max(LNUV,na.rm = T)) %>%
mutate(se=(Std_Dev^2)/sqrt(N)))## # A tibble: 6 x 9
## # Groups: GROUP [2]
## GROUP TIME Obs N Mean Std_Dev Minimum Maximum se
## <fct> <fct> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 control LNUV0 15 15 7.78 0.818 6.64 9.40 0.173
## 2 control LNUV4 15 15 7.82 0.784 6.37 9.41 0.159
## 3 control LNUV8 15 12 7.86 0.814 6.51 9.44 0.191
## 4 treatment LNUV0 19 19 7.59 0.449 6.96 8.39 0.0463
## 5 treatment LNUV4 19 19 7.42 0.503 6.53 8.27 0.0581
## 6 treatment LNUV8 19 15 7.40 0.455 6.88 8.52 0.0534ggplot(des, aes(x=TIME, y=Mean, group=GROUP, col=GROUP)) +
labs(x="Time (Months)", y="Log Urine Volume") +
geom_point(shape = 19, size = 2) +
geom_line() +
geom_errorbar(aes(ymin=Mean-se, ymax=Mean+se), width=0.25) +
theme_bw()
Pearson Correlation Coefficients
dat %>% dplyr::select(LNUV0, LNUV4, LNUV8) %>% cor(use="pairwise.complete.obs") # looks like compount symetry or Teoplitz## LNUV0 LNUV4 LNUV8
## LNUV0 1.0000000 0.9207956 0.9068861
## LNUV4 0.9207956 1.0000000 0.9171384
## LNUV8 0.9068861 0.9171384 1.0000000dat %>% dplyr::select(LNUV0, LNUV4, LNUV8) %>% ggcorr(label = TRUE, label_round = 4)
The Mixed Procedure
dt$GROUP<-fct_rev(dt$GROUP)
dt$TIME<-fct_rev(dt$TIME)
set_default_contrasts()
#fit.toep <- lme(LNUV ~ GROUP + TIME + GROUP*TIME, random = ~ 1 | PATIENT, data=dt, na.action=na.omit)
fit.toep <- lme(LNUV ~ GROUP + TIME + GROUP*TIME, random = ~ 1 | PATIENT, data=dt, na.action=na.omit,
corr = corARMA(form = ~ 1 | PATIENT, p = 0, q = 2))
summary(fit.toep)## Linear mixed-effects model fit by REML
## Data: dt
## AIC BIC logLik
## 99.00808 123.8944 -39.50404
##
## Random effects:
## Formula: ~1 | PATIENT
## (Intercept) Residual
## StdDev: 0.6283438 0.1851994
##
## Correlation Structure: ARMA(0,2)
## Formula: ~1 | PATIENT
## Parameter estimate(s):
## Theta1 Theta2
## 0.08069614 -0.10935064
## Fixed effects: LNUV ~ GROUP + TIME + GROUP * TIME
## Value Std.Error DF t-value p-value
## (Intercept) 7.453892 0.15271318 57 48.80975 0.0000
## GROUPcontrol 0.273125 0.22978980 32 1.18858 0.2433
## TIMELNUV4 -0.029905 0.06396836 57 -0.46750 0.6419
## TIMELNUV0 0.139410 0.06880678 57 2.02610 0.0474
## GROUPcontrol:TIMELNUV4 0.118654 0.09600415 57 1.23592 0.2216
## GROUPcontrol:TIMELNUV0 -0.084335 0.10330998 57 -0.81633 0.4177
## Correlation:
## (Intr) GROUPc TIMELNUV4 TIMELNUV0 GROUP:TIMELNUV4
## GROUPcontrol -0.665
## TIMELNUV4 -0.247 0.164
## TIMELNUV0 -0.260 0.173 0.621
## GROUPcontrol:TIMELNUV4 0.165 -0.245 -0.666 -0.414
## GROUPcontrol:TIMELNUV0 0.173 -0.259 -0.414 -0.666 0.619
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.81248221 -0.58170217 0.06693973 0.49265676 2.30805919
##
## Number of Observations: 95
## Number of Groups: 34Type III Analysis of Effects
set_sum_contrasts() # use sum coding, necessary to make type III LR tests valid
Anova(fit.toep,type="III")## Analysis of Deviance Table (Type III tests)
##
## Response: LNUV
## Chisq Df Pr(>Chisq)
## (Intercept) 2382.3914 1 < 2e-16 ***
## GROUP 1.4127 1 0.23460
## TIME 8.9624 2 0.01132 *
## GROUP:TIME 5.5852 2 0.06126 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Random coefficients (intercept and slope) model
Go to SAS output, pp 154-156, 157-161
head(dat)## SUBJECT GENDER AGE08 AGE10 AGE12 AGE14
## 1 1 1 21.0 20.0 21.5 23.0
## 2 2 1 21.0 21.5 24.0 25.5
## 3 3 1 20.5 24.0 24.5 26.0
## 4 4 1 23.5 24.5 25.0 26.5
## 5 5 1 21.5 23.0 22.5 23.5
## 6 6 1 20.0 21.0 21.0 22.5dt<-dat %>% gather(key = AGE, value = GROWTH, AGE08:AGE14) %>%
mutate(AGE=as.numeric(substr(AGE,4,5)))The MEANS procedure
(des<-dt %>%
group_by(GENDER, AGE) %>% summarize(N=sum(!is.na(GROWTH)), Mean=mean(GROWTH,na.rm = T), Std_Dev=sd(GROWTH,na.rm = T), Minimum=min(GROWTH,na.rm = T), Maximum=max(GROWTH,na.rm = T)) %>%
mutate(se=(Std_Dev^2)/sqrt(N)))## # A tibble: 8 x 8
## # Groups: GENDER [2]
## GENDER AGE N Mean Std_Dev Minimum Maximum se
## <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 8 16 22.9 2.45 17 27.5 1.50
## 2 0 10 16 23.8 2.14 20.5 28 1.14
## 3 0 12 16 25.7 2.65 22.5 31 1.76
## 4 0 14 16 27.5 2.09 25 31.5 1.09
## 5 1 8 11 21.2 2.12 16.5 24.5 1.36
## 6 1 10 11 22.2 1.90 19 25 1.09
## 7 1 12 11 23.1 2.36 19 28 1.69
## 8 1 14 11 24.1 2.44 19.5 28 1.79ggplot(des, aes(x=AGE, y=Mean, group=GENDER, col=factor(GENDER))) +
labs(x="AGE", y="GROWTH") +
geom_point(shape = 19, size = 2) +
geom_line() +
geom_errorbar(aes(ymin=Mean-se, ymax=Mean+se), width=0.25) +
theme_bw()
The Mixed Procedure
set_default_contrasts()
fit <- lme(GROWTH ~ AGE + GENDER + GENDER*AGE, random = ~ AGE | SUBJECT, data=dt)
summary(fit)## Linear mixed-effects model fit by REML
## Data: dt
## AIC BIC logLik
## 448.5817 469.7368 -216.2908
##
## Random effects:
## Formula: ~AGE | SUBJECT
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 2.4055009 (Intr)
## AGE 0.1803455 -0.668
## Residual 1.3100396
##
## Fixed effects: GROWTH ~ AGE + GENDER + GENDER * AGE
## Value Std.Error DF t-value p-value
## (Intercept) 16.340625 1.0185320 79 16.043311 0.0000
## AGE 0.784375 0.0859995 79 9.120691 0.0000
## GENDER 1.032102 1.5957329 25 0.646789 0.5237
## AGE:GENDER -0.304830 0.1347353 79 -2.262432 0.0264
## Correlation:
## (Intr) AGE GENDER
## AGE -0.880
## GENDER -0.638 0.562
## AGE:GENDER 0.562 -0.638 -0.880
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.168077730 -0.385939008 0.007104087 0.445154544 3.849463576
##
## Number of Observations: 108
## Number of Groups: 27Type III Analysis of Effects
set_sum_contrasts() # use sum coding, necessary to make type III LR tests valid
Anova(fit,type="III")## Analysis of Deviance Table (Type III tests)
##
## Response: GROWTH
## Chisq Df Pr(>Chisq)
## (Intercept) 257.3878 1 < 2e-16 ***
## AGE 83.1870 1 < 2e-16 ***
## GENDER 0.4183 1 0.51777
## AGE:GENDER 5.1186 1 0.02367 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Repeated measures regression model (with time-varying covariates)
Go to SAS output, pp 162-167, 168-176
head(dat)## # A tibble: 6 x 16
## IDNO AGE SEX GROUP Q13_1 CESD_1 Q13_2 CESD_2 Q13_3 CESD_3 Q13_4
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 10001 64 0 0 1 8 3 0 2 3 NA
## 2 10003 71 0 0 3 9 3 10 2 7 3
## 3 10005 82 0 0 4 12 NA NA NA NA NA
## 4 10007 70 1 0 2 3 3 7 3 2 2
## 5 10010 70 0 0 5 1 5 1 5 0 4
## 6 10012 76 0 0 2 16 NA NA 3 8 NA
## # … with 5 more variables: CESD_4 <dbl>, Q13_5 <dbl>, CESD_5 <dbl>,
## # AGEG <fct>, AGEGRP <dbl>dt <- dat %>%
gather(-IDNO,-AGE,-SEX,-GROUP,-AGEG,-AGEGRP,key = test,value = score) %>%
separate(test, into = c("test", "TIME")) %>%
spread(key = "test", value = "score") %>%
rename(HEALTH = Q13) %>%
mutate(HLTHGRP=ifelse(is.na(HEALTH),NA,ifelse(HEALTH<3,1,0)))The CORR Procedure
dat %>% dplyr::select(CESD_1,CESD_2,CESD_3,CESD_4,CESD_5) %>% cor(use="pairwise.complete.obs") # looks like compount symetry or Teoplitz## CESD_1 CESD_2 CESD_3 CESD_4 CESD_5
## CESD_1 1.0000000 0.5652716 0.5255772 0.4357750 0.5129300
## CESD_2 0.5652716 1.0000000 0.6294628 0.5152046 0.5270910
## CESD_3 0.5255772 0.6294628 1.0000000 0.5724844 0.5948377
## CESD_4 0.4357750 0.5152046 0.5724844 1.0000000 0.5193598
## CESD_5 0.5129300 0.5270910 0.5948377 0.5193598 1.0000000dat %>% dplyr::select(CESD_1,CESD_2,CESD_3,CESD_4,CESD_5) %>% ggcorr(label = TRUE, label_round = 4)
The MEANS procedure
(des_sex<-dt %>%
group_by(SEX, TIME) %>% summarize(Obs=n(), N=sum(!is.na(CESD)), Mean=mean(CESD,na.rm = T), Std_Dev=sd(CESD,na.rm = T), Minimum=min(CESD,na.rm = T), Maximum=max(CESD,na.rm = T)) %>%
mutate(se=(Std_Dev^2)/sqrt(N)))## # A tibble: 10 x 9
## # Groups: SEX [2]
## SEX TIME Obs N Mean Std_Dev Minimum Maximum se
## <dbl> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 1 490 487 7.95 7.29 0 42 2.41
## 2 0 2 490 365 8.40 7.12 0 38 2.66
## 3 0 3 490 336 7.44 7.16 0 36 2.80
## 4 0 4 490 231 7.83 8.08 0 43 4.30
## 5 0 5 490 195 7.78 6.62 0 38 3.14
## 6 1 1 408 402 6.14 6.15 0 40 1.89
## 7 1 2 408 339 6.18 6.73 0 34 2.46
## 8 1 3 408 303 5.56 5.87 0 31 1.98
## 9 1 4 408 237 5.45 6.30 0 37 2.58
## 10 1 5 408 191 5.75 6.12 0 44 2.71(des_age<-dt %>%
group_by(AGEGRP, TIME) %>% summarize(Obs=n(), N=sum(!is.na(CESD)), Mean=mean(CESD,na.rm = T), Std_Dev=sd(CESD,na.rm = T), Minimum=min(CESD,na.rm = T), Maximum=max(CESD,na.rm = T)) %>%
mutate(se=(Std_Dev^2)/sqrt(N)))## # A tibble: 25 x 9
## # Groups: AGEGRP [5]
## AGEGRP TIME Obs N Mean Std_Dev Minimum Maximum se
## <dbl> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 264 263 6.57 6.70 0 37 2.77
## 2 1 2 264 211 7.44 7.51 0 38 3.89
## 3 1 3 264 195 6.48 6.71 0 32 3.22
## 4 1 4 264 153 6.90 7.24 0 37 4.24
## 5 1 5 264 130 6.59 6.88 0 38 4.15
## 6 2 1 322 319 7.04 7.14 0 42 2.85
## 7 2 2 322 261 7.51 7.08 0 33 3.11
## 8 2 3 322 242 6.53 6.93 0 36 3.09
## 9 2 4 322 184 6.54 7.80 0 43 4.48
## 10 2 5 322 155 7.30 7.01 0 44 3.95
## # … with 15 more rows(des_health<-dt %>%
group_by(HLTHGRP, TIME) %>% summarize(Obs=n(), N=sum(!is.na(CESD)), Mean=mean(CESD,na.rm = T), Std_Dev=sd(CESD,na.rm = T), Minimum=min(CESD,na.rm = T), Maximum=max(CESD,na.rm = T)) %>%
mutate(se=(Std_Dev^2)/sqrt(N)))## # A tibble: 15 x 9
## # Groups: HLTHGRP [3]
## HLTHGRP TIME Obs N Mean Std_Dev Minimum Maximum se
## <dbl> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0 1 746 739 6.41 6.16 0 34 1.40
## 2 0 2 581 577 6.60 6.23 0 35 1.62
## 3 0 3 575 573 5.95 6.10 0 36 1.55
## 4 0 4 414 411 5.90 6.70 0 43 2.22
## 5 0 5 344 343 6.41 6.06 0 44 1.98
## 6 1 1 132 132 11.3 8.87 0 42 6.85
## 7 1 2 112 111 10.9 9.28 0 38 8.17
## 8 1 3 60 60 11.7 8.83 0 32 10.1
## 9 1 4 54 54 12.1 9.48 0 37 12.2
## 10 1 5 45 42 9.5 8.52 0 30 11.2
## 11 NA 1 20 18 5.94 6.01 0 21 8.53
## 12 NA 2 205 16 8.75 8.25 0 22 17.0
## 13 NA 3 263 6 12.2 7.70 6 27 24.2
## 14 NA 4 430 3 7 3 4 10 5.20
## 15 NA 5 509 1 18 NA 18 18 NAThe Mixed Procedure
dt$TIME<-fct_rev(dt$TIME)
set_default_contrasts()
fit1 <- lme(CESD ~ AGE + SEX + HLTHGRP + TIME + AGE*TIME + SEX*TIME + HLTHGRP*TIME , random = ~ 1 | IDNO, data=dt, na.action=na.omit)
summary(fit1)## Linear mixed-effects model fit by REML
## Data: dt
## AIC BIC logLik
## 19299.21 19431.51 -9627.604
##
## Random effects:
## Formula: ~1 | IDNO
## (Intercept) Residual
## StdDev: 4.830165 4.614
##
## Fixed effects: CESD ~ AGE + SEX + HLTHGRP + TIME + AGE * TIME + SEX * TIME + HLTHGRP * TIME
## Value Std.Error DF t-value p-value
## (Intercept) 2.971781 4.802607 2133 0.618785 0.5361
## AGE 0.069914 0.065942 889 1.060236 0.2893
## SEX -2.029810 0.604267 889 -3.359124 0.0008
## HLTHGRP 1.745612 0.836775 2133 2.086118 0.0371
## TIME4 1.849726 5.172271 2133 0.357624 0.7207
## TIME3 1.692825 4.922739 2133 0.343879 0.7310
## TIME2 6.209829 4.934854 2133 1.258361 0.2084
## TIME1 -3.921048 4.777841 2133 -0.820674 0.4119
## AGE:TIME4 -0.029518 0.071161 2133 -0.414803 0.6783
## AGE:TIME3 -0.031908 0.067593 2133 -0.472066 0.6369
## AGE:TIME2 -0.085356 0.067810 2133 -1.258749 0.2083
## AGE:TIME1 0.045920 0.065615 2133 0.699830 0.4841
## SEX:TIME4 -0.371479 0.645365 2133 -0.575610 0.5649
## SEX:TIME3 0.075114 0.616741 2133 0.121793 0.9031
## SEX:TIME2 -0.210267 0.614191 2133 -0.342348 0.7321
## SEX:TIME1 0.367271 0.599012 2133 0.613127 0.5399
## HLTHGRP:TIME4 2.648663 1.072792 2133 2.468945 0.0136
## HLTHGRP:TIME3 0.937435 1.072002 2133 0.874472 0.3820
## HLTHGRP:TIME2 1.043784 0.981748 2133 1.063189 0.2878
## HLTHGRP:TIME1 0.964092 0.968151 2133 0.995808 0.3195
## Correlation:
## (Intr) AGE SEX HLTHGRP TIME4 TIME3 TIME2 TIME1
## AGE -0.996
## SEX -0.120 0.059
## HLTHGRP -0.051 0.031 0.015
## TIME4 -0.605 0.603 0.077 0.047
## TIME3 -0.669 0.668 0.079 0.047 0.587
## TIME2 -0.691 0.690 0.084 0.049 0.588 0.650
## TIME1 -0.744 0.743 0.088 0.050 0.608 0.672 0.693
## AGE:TIME4 0.602 -0.606 -0.040 -0.030 -0.996 -0.585 -0.585 -0.606
## AGE:TIME3 0.668 -0.673 -0.038 -0.029 -0.586 -0.996 -0.649 -0.671
## AGE:TIME2 0.690 -0.695 -0.043 -0.030 -0.586 -0.649 -0.996 -0.691
## AGE:TIME1 0.743 -0.748 -0.043 -0.031 -0.607 -0.671 -0.692 -0.996
## SEX:TIME4 0.077 -0.040 -0.594 -0.013 -0.114 -0.074 -0.074 -0.077
## SEX:TIME3 0.079 -0.038 -0.652 -0.011 -0.074 -0.127 -0.078 -0.079
## SEX:TIME2 0.085 -0.043 -0.674 -0.017 -0.075 -0.078 -0.117 -0.085
## SEX:TIME1 0.089 -0.044 -0.715 -0.013 -0.077 -0.079 -0.084 -0.119
## HLTHGRP:TIME4 0.038 -0.025 -0.009 -0.715 -0.034 -0.038 -0.036 -0.039
## HLTHGRP:TIME3 0.038 -0.024 -0.008 -0.735 -0.036 -0.055 -0.034 -0.040
## HLTHGRP:TIME2 0.047 -0.032 -0.010 -0.820 -0.037 -0.044 -0.050 -0.050
## HLTHGRP:TIME1 0.044 -0.029 -0.006 -0.837 -0.039 -0.045 -0.041 -0.060
## AGE:TIME4 AGE:TIME3 AGE:TIME2 AGE:TIME1 SEX:TIME4 SEX:TIME3
## AGE
## SEX
## HLTHGRP
## TIME4
## TIME3
## TIME2
## TIME1
## AGE:TIME4
## AGE:TIME3 0.589
## AGE:TIME2 0.588 0.653
## AGE:TIME1 0.609 0.676 0.696
## SEX:TIME4 0.051 0.039 0.039 0.041
## SEX:TIME3 0.038 0.066 0.038 0.038 0.581
## SEX:TIME2 0.039 0.039 0.057 0.043 0.583 0.638
## SEX:TIME1 0.040 0.038 0.042 0.057 0.599 0.657
## HLTHGRP:TIME4 0.010 0.025 0.022 0.025 0.025 0.011
## HLTHGRP:TIME3 0.023 0.030 0.020 0.025 0.011 0.048
## HLTHGRP:TIME2 0.023 0.029 0.026 0.034 0.012 0.011
## HLTHGRP:TIME1 0.025 0.029 0.025 0.035 0.010 0.008
## SEX:TIME2 SEX:TIME1 HLTHGRP:TIME4 HLTHGRP:TIME3
## AGE
## SEX
## HLTHGRP
## TIME4
## TIME3
## TIME2
## TIME1
## AGE:TIME4
## AGE:TIME3
## AGE:TIME2
## AGE:TIME1
## SEX:TIME4
## SEX:TIME3
## SEX:TIME2
## SEX:TIME1 0.678
## HLTHGRP:TIME4 0.010 0.012
## HLTHGRP:TIME3 0.009 0.009 0.581
## HLTHGRP:TIME2 -0.002 0.014 0.624 0.654
## HLTHGRP:TIME1 0.010 0.031 0.635 0.661
## HLTHGRP:TIME2
## AGE
## SEX
## HLTHGRP
## TIME4
## TIME3
## TIME2
## TIME1
## AGE:TIME4
## AGE:TIME3
## AGE:TIME2
## AGE:TIME1
## SEX:TIME4
## SEX:TIME3
## SEX:TIME2
## SEX:TIME1
## HLTHGRP:TIME4
## HLTHGRP:TIME3
## HLTHGRP:TIME2
## HLTHGRP:TIME1 0.728
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.5348966 -0.5419375 -0.1354727 0.4064318 6.9225994
##
## Number of Observations: 3042
## Number of Groups: 892Type III Analysis of Effects
set_sum_contrasts() # use sum coding, necessary to make type III LR tests valid
Anova(fit1,type="III")## Analysis of Deviance Table (Type III tests)
##
## Response: CESD
## Chisq Df Pr(>Chisq)
## (Intercept) 0.3829 1 0.5360580
## AGE 1.1241 1 0.2890373
## SEX 11.2837 1 0.0007819 ***
## HLTHGRP 4.3519 1 0.0369680 *
## TIME 7.2964 4 0.1210286
## AGE:TIME 6.6944 4 0.1529462
## SEX:TIME 2.3059 4 0.6797011
## HLTHGRP:TIME 6.8161 4 0.1459337
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Test covariance structures
head(dt)## PREPARATION SUBJECT TIME CAR TIME.N
## 1 1 71 CAR_00 116 0
## 2 1 73 CAR_00 146 0
## 3 1 80 CAR_00 200 0
## 4 1 83 CAR_00 180 0
## 5 1 90 CAR_00 142 0
## 6 1 92 CAR_00 106 0Pearson Correlation Coefficients
dat %>% dplyr::select(CAR_00, CAR_06, CAR_12) %>% cor(use="pairwise.complete.obs") # looks like compount symetry or Teoplitz## CAR_00 CAR_06 CAR_12
## CAR_00 1.0000000 0.6509899 0.7110154
## CAR_06 0.6509899 1.0000000 0.8633590
## CAR_12 0.7110154 0.8633590 1.0000000dat %>% dplyr::select(CAR_00, CAR_06, CAR_12) %>% ggcorr(label = TRUE, label_round = 4)
The Mixed Procedure
set_default_contrasts()
fit.sym <- lme(CAR ~ TIME.N + PREPARATION + PREPARATION*TIME.N, random = ~ 1 | SUBJECT, data=dt, na.action=na.omit,
corr = corCompSymm(form = ~ 1 | SUBJECT))
fit.ar <- lme(CAR ~ TIME.N + PREPARATION + PREPARATION*TIME.N, random = ~ 1 | SUBJECT, data=dt, na.action=na.omit,
corr = corAR1(form = ~ 1 | SUBJECT))
fit.toep <- lme(CAR ~ TIME.N + PREPARATION + PREPARATION*TIME.N, random = ~ 1 | SUBJECT, data=dt, na.action=na.omit,
corr = corARMA(form = ~ 1 | SUBJECT, p = 0, q = 2))
fit.un <- lme(CAR ~ TIME.N + PREPARATION + PREPARATION*TIME.N, random = ~ 1 | SUBJECT, data=dt, na.action=na.omit)
summary(fit.sym)## Linear mixed-effects model fit by REML
## Data: dt
## AIC BIC logLik
## 552.2793 568.5393 -267.1397
##
## Random effects:
## Formula: ~1 | SUBJECT
## (Intercept) Residual
## StdDev: 62.15739 52.88666
##
## Correlation Structure: Compound symmetry
## Formula: ~1 | SUBJECT
## Parameter estimate(s):
## Rho
## 0
## Fixed effects: CAR ~ TIME.N + PREPARATION + PREPARATION * TIME.N
## Value Std.Error DF t-value p-value
## (Intercept) 158.38889 32.13091 31 4.929486 0.0000
## TIME.N 10.36111 2.54451 31 4.071946 0.0003
## PREPARATION2 -21.66667 45.43997 14 -0.476820 0.6408
## PREPARATION3 22.54444 47.65785 14 0.473048 0.6435
## TIME.N:PREPARATION2 -4.83333 3.59848 31 -1.343159 0.1890
## TIME.N:PREPARATION3 6.03889 3.77412 31 1.600079 0.1197
## Correlation:
## (Intr) TIME.N PREPARATION2 PREPARATION3
## TIME.N -0.475
## PREPARATION2 -0.707 0.336
## PREPARATION3 -0.674 0.320 0.477
## TIME.N:PREPARATION2 0.336 -0.707 -0.475 -0.227
## TIME.N:PREPARATION3 0.320 -0.674 -0.227 -0.475
## TIME.N:PREPARATION2
## TIME.N
## PREPARATION2
## PREPARATION3
## TIME.N:PREPARATION2
## TIME.N:PREPARATION3 0.477
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.5156694 -0.5615791 -0.1263106 0.3838318 2.3883377
##
## Number of Observations: 51
## Number of Groups: 17summary(fit.ar)## Linear mixed-effects model fit by REML
## Data: dt
## AIC BIC logLik
## 552.2787 568.5386 -267.1393
##
## Random effects:
## Formula: ~1 | SUBJECT
## (Intercept) Residual
## StdDev: 62.44697 52.62844
##
## Correlation Structure: AR(1)
## Formula: ~1 | SUBJECT
## Parameter estimate(s):
## Phi
## -0.01191066
## Fixed effects: CAR ~ TIME.N + PREPARATION + PREPARATION * TIME.N
## Value Std.Error DF t-value p-value
## (Intercept) 158.46842 32.12722 31 4.932528 0.0000
## TIME.N 10.36111 2.53191 31 4.092216 0.0003
## PREPARATION2 -21.61658 45.43476 14 -0.475772 0.6416
## PREPARATION3 22.60517 47.65237 14 0.474377 0.6425
## TIME.N:PREPARATION2 -4.83333 3.58066 31 -1.349845 0.1868
## TIME.N:PREPARATION3 6.03889 3.75543 31 1.608044 0.1180
## Correlation:
## (Intr) TIME.N PREPARATION2 PREPARATION3
## TIME.N -0.473
## PREPARATION2 -0.707 0.334
## PREPARATION3 -0.674 0.319 0.477
## TIME.N:PREPARATION2 0.334 -0.707 -0.473 -0.225
## TIME.N:PREPARATION3 0.319 -0.674 -0.225 -0.473
## TIME.N:PREPARATION2
## TIME.N
## PREPARATION2
## PREPARATION3
## TIME.N:PREPARATION2
## TIME.N:PREPARATION3 0.477
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.5361194 -0.5643613 -0.1250449 0.3855594 2.3838833
##
## Number of Observations: 51
## Number of Groups: 17summary(fit.toep)## Linear mixed-effects model fit by REML
## Data: dt
## AIC BIC logLik
## 554.2787 572.3453 -267.1393
##
## Random effects:
## Formula: ~1 | SUBJECT
## (Intercept) Residual
## StdDev: 62.30788 52.79374
##
## Correlation Structure: ARMA(0,2)
## Formula: ~1 | SUBJECT
## Parameter estimate(s):
## Theta1 Theta2
## -0.005546214 0.006419251
## Fixed effects: CAR ~ TIME.N + PREPARATION + PREPARATION * TIME.N
## Value Std.Error DF t-value p-value
## (Intercept) 158.46857 32.12737 31 4.932510 0.0000
## TIME.N 10.36111 2.53187 31 4.092268 0.0003
## PREPARATION2 -21.61648 45.43496 14 -0.475768 0.6416
## PREPARATION3 22.60528 47.65259 14 0.474377 0.6425
## TIME.N:PREPARATION2 -4.83333 3.58061 31 -1.349862 0.1868
## TIME.N:PREPARATION3 6.03889 3.75538 31 1.608064 0.1180
## Correlation:
## (Intr) TIME.N PREPARATION2 PREPARATION3
## TIME.N -0.473
## PREPARATION2 -0.707 0.334
## PREPARATION3 -0.674 0.319 0.477
## TIME.N:PREPARATION2 0.334 -0.707 -0.473 -0.225
## TIME.N:PREPARATION3 0.319 -0.674 -0.225 -0.473
## TIME.N:PREPARATION2
## TIME.N
## PREPARATION2
## PREPARATION3
## TIME.N:PREPARATION2
## TIME.N:PREPARATION3 0.477
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.5260907 -0.5647578 -0.1278987 0.3840093 2.3824212
##
## Number of Observations: 51
## Number of Groups: 17summary(fit.un)## Linear mixed-effects model fit by REML
## Data: dt
## AIC BIC logLik
## 550.2793 564.7326 -267.1397
##
## Random effects:
## Formula: ~1 | SUBJECT
## (Intercept) Residual
## StdDev: 62.15739 52.88666
##
## Fixed effects: CAR ~ TIME.N + PREPARATION + PREPARATION * TIME.N
## Value Std.Error DF t-value p-value
## (Intercept) 158.38889 32.13091 31 4.929486 0.0000
## TIME.N 10.36111 2.54451 31 4.071946 0.0003
## PREPARATION2 -21.66667 45.43997 14 -0.476820 0.6408
## PREPARATION3 22.54444 47.65785 14 0.473048 0.6435
## TIME.N:PREPARATION2 -4.83333 3.59848 31 -1.343159 0.1890
## TIME.N:PREPARATION3 6.03889 3.77412 31 1.600079 0.1197
## Correlation:
## (Intr) TIME.N PREPARATION2 PREPARATION3
## TIME.N -0.475
## PREPARATION2 -0.707 0.336
## PREPARATION3 -0.674 0.320 0.477
## TIME.N:PREPARATION2 0.336 -0.707 -0.475 -0.227
## TIME.N:PREPARATION3 0.320 -0.674 -0.227 -0.475
## TIME.N:PREPARATION2
## TIME.N
## PREPARATION2
## PREPARATION3
## TIME.N:PREPARATION2
## TIME.N:PREPARATION3 0.477
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.5156694 -0.5615791 -0.1263106 0.3838318 2.3883377
##
## Number of Observations: 51
## Number of Groups: 17AIC(fit.sym, fit.ar, fit.toep, fit.un)## df AIC
## fit.sym 9 552.2793
## fit.ar 9 552.2787
## fit.toep 10 554.2787
## fit.un 8 550.2793Unstructured resulting in smallest AIC is considered best.