[1] "C:/Users/jsteibel/OneDrive/Documents/job/500 ANS/2025/14"Week 14
Repeated Measures
Setup Code
#================================================================#
# Setup Options
#================================================================#
# good habit to run is the sessionInfo() function to print all relevant info
sessionInfo()R version 4.4.2 (2024-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 10 x64 (build 19045)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.utf8
[2] LC_CTYPE=English_United States.utf8
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.utf8
time zone: America/Chicago
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] htmlwidgets_1.6.4 compiler_4.4.2 fastmap_1.2.0 cli_3.6.3
[5] tools_4.4.2 htmltools_0.5.8.1 rstudioapi_0.17.1 yaml_2.3.10
[9] rmarkdown_2.29 knitr_1.49 jsonlite_1.8.9 xfun_0.49
[13] digest_0.6.37 rlang_1.1.4 evaluate_1.0.1 # remove all objects if restarting script
rm(list=ls())
knitr::opts_knit$set(root.dir = "C:/Users/jsteibel/OneDrive/Documents/job/500 ANS/2025/14/")
getwd()[1] "C:/Users/jsteibel/OneDrive/Documents/job/500 ANS/2025/14"# list all objects in your environment
ls()character(0)# set tibble width for printing - print all columns
options(tibble.width = Inf)
# remove scientific notation (by default will print)
options(scipen=999)
#================================================================#
# Packages
#================================================================#
# Here we must load the packages for the current environment to make them
# accessible.
# load libraries
library(tidyverse)
library(tidymodels)
library(emmeans)
library(car)
library(DT)
library(gridExtra)
library(lme4)
library(lmerTest)
library(SASmixed)
library(mmrm)Data input 1
First dataset: read weights from 30 rats over several weeks (days of age is the time variable)
workdir="C:/Users/jsteibel/OneDrive/Documents/job/500 ANS/2025/14/"
rats<-read.csv(paste0(workdir,"rats_rm.csv"))[,-1]
rats indiv d_8 d_15 d_22 d_29 d_36
1 1 151 199 246 283 320
2 2 145 199 249 293 354
3 3 147 214 263 312 328
4 4 155 200 237 272 297
5 5 135 188 230 280 323
6 6 159 210 252 298 331
7 7 141 189 231 275 305
8 8 159 201 248 297 338
9 9 177 236 285 350 376
10 10 134 182 220 260 296
11 11 160 208 261 313 352
12 12 143 188 220 273 314
13 13 154 200 244 289 325
14 14 171 221 270 326 358
15 15 163 216 242 281 312
16 16 160 207 248 288 324
17 17 142 187 234 280 316
18 18 156 203 243 283 317
19 19 157 212 259 307 336
20 20 152 203 246 286 321
21 21 154 205 253 298 334
22 22 139 190 225 267 302
23 23 146 191 229 272 302
24 24 157 211 250 285 323
25 25 132 185 237 286 331
26 26 160 207 257 303 345
27 27 169 216 261 295 333
28 28 157 205 248 289 316
29 29 137 180 219 258 291
30 30 153 200 244 286 324glimpse(rats)Rows: 30
Columns: 6
$ indiv <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 1…
$ d_8 <int> 151, 145, 147, 155, 135, 159, 141, 159, 177, 134, 160, 143, 154,…
$ d_15 <int> 199, 199, 214, 200, 188, 210, 189, 201, 236, 182, 208, 188, 200,…
$ d_22 <int> 246, 249, 263, 237, 230, 252, 231, 248, 285, 220, 261, 220, 244,…
$ d_29 <int> 283, 293, 312, 272, 280, 298, 275, 297, 350, 260, 313, 273, 289,…
$ d_36 <int> 320, 354, 328, 297, 323, 331, 305, 338, 376, 296, 352, 314, 325,…long_rats<-pivot_longer(data = rats,cols = d_8:d_36,names_to="days",names_prefix = "d_",values_to= "weight")
long_rats<-mutate(long_rats,indiv=as.factor(indiv),days=as.numeric(days))
long_rats# A tibble: 150 × 3
indiv days weight
<fct> <dbl> <int>
1 1 8 151
2 1 15 199
3 1 22 246
4 1 29 283
5 1 36 320
6 2 8 145
7 2 15 199
8 2 22 249
9 2 29 293
10 2 36 354
# ℹ 140 more rowsRepresent the data
ggplot(long_rats,aes(x=days,y=weight,color=indiv))+geom_point()+geom_line()+geom_smooth(aes(color=NULL))`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Start fitting models
Saturated mean and variance model
long_rats<-mutate(long_rats,period=as.factor(days))
long_rats# A tibble: 150 × 4
indiv days weight period
<fct> <dbl> <int> <fct>
1 1 8 151 8
2 1 15 199 15
3 1 22 246 22
4 1 29 283 29
5 1 36 320 36
6 2 8 145 8
7 2 15 199 15
8 2 22 249 22
9 2 29 293 29
10 2 36 354 36
# ℹ 140 more rowsfit_full<-mmrm(formula=weight~period-1+us(period|indiv),data = long_rats)
fit_fullmmrm fit
Formula: weight ~ period - 1 + us(period | indiv)
Data: long_rats (used 150 observations from 30 subjects with maximum 5
timepoints)
Covariance: unstructured (15 variance parameters)
Inference: REML
Deviance: 978.7812
Coefficients:
period8 period15 period22 period29 period36
152.1667 201.7667 245.0333 289.5000 324.8000
Model Inference Optimization:
Converged with code 0 and message: convergence: rel_reduction_of_f <= factr*epsmchsummary(fit_full)mmrm fit
Formula: weight ~ period - 1 + us(period | indiv)
Data: long_rats (used 150 observations from 30 subjects with maximum 5
timepoints)
Covariance: unstructured (15 variance parameters)
Method: Satterthwaite
Vcov Method: Asymptotic
Inference: REML
Model selection criteria:
AIC BIC logLik deviance
1008.8 1029.8 -489.4 978.8
Coefficients:
Estimate Std. Error df t value Pr(>|t|)
period8 152.167 2.038 29.027 74.68 <0.0000000000000002 ***
period15 201.767 2.313 29.030 87.24 <0.0000000000000002 ***
period22 245.033 2.805 29.039 87.36 <0.0000000000000002 ***
period29 289.500 3.496 29.043 82.81 <0.0000000000000002 ***
period36 324.800 3.551 29.042 91.46 <0.0000000000000002 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Covariance estimate:
8 15 22 29 36
8 124.5663 129.6434 141.3882 155.1270 133.5079
15 129.6434 160.4727 179.4888 204.6157 172.7079
22 141.3882 179.4888 236.0054 278.7911 256.4632
29 155.1270 204.6157 278.7911 366.6797 342.5501
36 133.5079 172.7079 256.4632 342.5501 378.3175VarCorr(fit_full) 8 15 22 29 36
8 124.5663 129.6434 141.3882 155.1270 133.5079
15 129.6434 160.4727 179.4888 204.6157 172.7079
22 141.3882 179.4888 236.0054 278.7911 256.4632
29 155.1270 204.6157 278.7911 366.6797 342.5501
36 133.5079 172.7079 256.4632 342.5501 378.3175cov2cor(VarCorr(fit_full)) 8 15 22 29 36
8 1.0000000 0.9169581 0.8246165 0.7258446 0.6150049
15 0.9169581 1.0000000 0.9223075 0.8435189 0.7009441
22 0.8246165 0.9223075 1.0000000 0.9477078 0.8582935
29 0.7258446 0.8435189 0.9477078 1.0000000 0.9197132
36 0.6150049 0.7009441 0.8582935 0.9197132 1.0000000#cov? #cor?
cr<-VarCorr(fit_full)
cr<-cr%>%as.tibble()%>%mutate(ID=rownames(cr))Warning: `as.tibble()` was deprecated in tibble 2.0.0.
ℹ Please use `as_tibble()` instead.
ℹ The signature and semantics have changed, see `?as_tibble`.(cr)# A tibble: 5 × 6
`8` `15` `22` `29` `36` ID
<dbl> <dbl> <dbl> <dbl> <dbl> <chr>
1 125. 130. 141. 155. 134. 8
2 130. 160. 179. 205. 173. 15
3 141. 179. 236. 279. 256. 22
4 155. 205. 279. 367. 343. 29
5 134. 173. 256. 343. 378. 36 cr<-pivot_longer(cr,cols = 1:5,names_to="ID2",values_to= "cov")%>%mutate(ID=as.numeric(ID),ID2=as.numeric(ID2))
p1<-ggplot(cr,aes(ID,ID2,z=cov, ))+geom_contour()+stat_contour(geom="polygon",aes(fill=stat(level))) +
scale_fill_distiller(palette = "Spectral", direction = -1)
print(p1)Warning: `stat(level)` was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(level)` instead.
Saturated mean and variance cs
fit_cs<-mmrm(formula=weight~period-1+cs(period|indiv),data = long_rats)
summary(fit_cs)mmrm fit
Formula: weight ~ period - 1 + cs(period | indiv)
Data: long_rats (used 150 observations from 30 subjects with maximum 5
timepoints)
Covariance: compound symmetry (2 variance parameters)
Method: Satterthwaite
Vcov Method: Asymptotic
Inference: REML
Model selection criteria:
AIC BIC logLik deviance
1096.5 1099.4 -546.3 1092.5
Coefficients:
Estimate Std. Error df t value Pr(>|t|)
period8 152.167 2.906 41.640 52.35 <0.0000000000000002 ***
period15 201.767 2.906 41.640 69.42 <0.0000000000000002 ***
period22 245.033 2.906 41.640 84.31 <0.0000000000000002 ***
period29 289.500 2.906 41.640 99.61 <0.0000000000000002 ***
period36 324.800 2.906 41.640 111.75 <0.0000000000000002 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Covariance estimate:
8 15 22 29 36
8 253.4262 199.6432 199.6432 199.6432 199.6432
15 199.6432 253.4262 199.6432 199.6432 199.6432
22 199.6432 199.6432 253.4262 199.6432 199.6432
29 199.6432 199.6432 199.6432 253.4262 199.6432
36 199.6432 199.6432 199.6432 199.6432 253.4262cov2cor(VarCorr(fit_cs)) 8 15 22 29 36
8 1.0000000 0.7877765 0.7877765 0.7877765 0.7877765
15 0.7877765 1.0000000 0.7877765 0.7877765 0.7877765
22 0.7877765 0.7877765 1.0000000 0.7877765 0.7877765
29 0.7877765 0.7877765 0.7877765 1.0000000 0.7877765
36 0.7877765 0.7877765 0.7877765 0.7877765 1.0000000Saturated mean and variance csh
fit_csh<-mmrm(formula=weight~period-1+csh(period|indiv),data = long_rats)
summary(fit_csh)mmrm fit
Formula: weight ~ period - 1 + csh(period | indiv)
Data: long_rats (used 150 observations from 30 subjects with maximum 5
timepoints)
Covariance: heterogeneous compound symmetry (6 variance parameters)
Method: Satterthwaite
Vcov Method: Asymptotic
Inference: REML
Model selection criteria:
AIC BIC logLik deviance
1067.0 1075.5 -527.5 1055.0
Coefficients:
Estimate Std. Error df t value Pr(>|t|)
period8 152.167 2.120 28.493 71.77 <0.0000000000000002 ***
period15 201.767 2.287 29.507 88.22 <0.0000000000000002 ***
period22 245.033 2.698 30.050 90.81 <0.0000000000000002 ***
period29 289.500 3.427 29.683 84.49 <0.0000000000000002 ***
period36 324.800 3.686 28.546 88.11 <0.0000000000000002 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Covariance estimate:
8 15 22 29 36
8 134.8775 120.6886 142.3932 180.8293 194.5291
15 120.6886 156.9040 153.5806 195.0366 209.8127
22 142.3932 153.5806 218.4135 230.1117 247.5452
29 180.8293 195.0366 230.1117 352.2400 314.3650
36 194.5291 209.8127 247.5452 314.3650 407.6339cov2cor(VarCorr(fit_csh)) 8 15 22 29 36
8 1.0000000 0.8296208 0.8296208 0.8296208 0.8296208
15 0.8296208 1.0000000 0.8296208 0.8296208 0.8296208
22 0.8296208 0.8296208 1.0000000 0.8296208 0.8296208
29 0.8296208 0.8296208 0.8296208 1.0000000 0.8296208
36 0.8296208 0.8296208 0.8296208 0.8296208 1.0000000Saturated mean and variance ante1
fit_ante1<-mmrm(formula=weight~period-1+adh(period|indiv),data = long_rats)
summary(fit_ante1)mmrm fit
Formula: weight ~ period - 1 + adh(period | indiv)
Data: long_rats (used 150 observations from 30 subjects with maximum 5
timepoints)
Covariance: heterogeneous ante-dependence (9 variance parameters)
Method: Satterthwaite
Vcov Method: Asymptotic
Inference: REML
Model selection criteria:
AIC BIC logLik deviance
1006.6 1019.2 -494.3 988.6
Coefficients:
Estimate Std. Error df t value Pr(>|t|)
period8 152.167 2.038 29.008 74.67 <0.0000000000000002 ***
period15 201.767 2.313 29.011 87.22 <0.0000000000000002 ***
period22 245.033 2.805 29.015 87.34 <0.0000000000000002 ***
period29 289.500 3.497 29.017 82.79 <0.0000000000000002 ***
period36 324.800 3.552 29.018 91.45 <0.0000000000000002 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Covariance estimate:
8 15 22 29 36
8 124.6022 129.6968 145.0788 171.3822 160.1093
15 129.6968 160.5460 179.5866 212.1464 198.1922
22 145.0788 179.5866 236.1208 278.9304 260.5834
29 171.3822 212.1464 278.9304 366.8300 342.7013
36 160.1093 198.1922 260.5834 342.7013 378.4567cov2cor(VarCorr(fit_ante1)) 8 15 22 29 36
8 1.0000000 0.9169945 0.8458123 0.8016233 0.7373022
15 0.9169945 1.0000000 0.9223745 0.8741855 0.8040421
22 0.8458123 0.9223745 1.0000000 0.9477556 0.8717090
29 0.8016233 0.8741855 0.9477556 1.0000000 0.9197614
36 0.7373022 0.8040421 0.8717090 0.9197614 1.0000000Saturated mean and variance toeplitz
fit_toep<-mmrm(formula=weight~period-1+toeph(period|indiv),data = long_rats)Warning in mmrm(formula = weight ~ period - 1 + toeph(period | indiv), data =
long_rats): Divergence with optimizer L-BFGS-B due to problems: L-BFGS-B needs
finite values of 'fn'summary(fit_toep)mmrm fit
Formula: weight ~ period - 1 + toeph(period | indiv)
Data: long_rats (used 150 observations from 30 subjects with maximum 5
timepoints)
Covariance: heterogeneous Toeplitz (9 variance parameters)
Method: Satterthwaite
Vcov Method: Asymptotic
Inference: REML
Model selection criteria:
AIC BIC logLik deviance
999.7 1012.3 -490.9 981.7
Coefficients:
Estimate Std. Error df t value Pr(>|t|)
period8 152.167 2.080 32.410 73.15 <0.0000000000000002 ***
period15 201.767 2.317 31.550 87.09 <0.0000000000000002 ***
period22 245.033 2.756 30.160 88.92 <0.0000000000000002 ***
period29 289.500 3.417 31.860 84.73 <0.0000000000000002 ***
period36 324.800 3.475 33.540 93.46 <0.0000000000000002 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Covariance estimate:
8 15 22 29 36
8 129.8289 133.7940 143.1982 150.0235 129.8438
15 133.7940 161.0305 177.2309 197.7460 169.9296
22 143.1982 177.2309 227.8123 261.3821 239.2124
29 150.0235 197.7460 261.3821 350.2523 329.6243
36 129.8438 169.9296 239.2124 329.6243 362.2962cov2cor(VarCorr(fit_toep)) 8 15 22 29 36
8 1.0000000 0.9253303 0.8326513 0.7035309 0.5986920
15 0.9253303 1.0000000 0.9253303 0.8326513 0.7035309
22 0.8326513 0.9253303 1.0000000 0.9253303 0.8326513
29 0.7035309 0.8326513 0.9253303 1.0000000 0.9253303
36 0.5986920 0.7035309 0.8326513 0.9253303 1.0000000Saturated mean and variance autoregressive heterogeneous
fit_ar1h<-mmrm(formula=weight~period-1+ar1h(period|indiv),data = long_rats)
summary(fit_ar1h)mmrm fit
Formula: weight ~ period - 1 + ar1h(period | indiv)
Data: long_rats (used 150 observations from 30 subjects with maximum 5
timepoints)
Covariance: heterogeneous auto-regressive order one (6 variance parameters)
Method: Satterthwaite
Vcov Method: Asymptotic
Inference: REML
Model selection criteria:
AIC BIC logLik deviance
1002.4 1010.8 -495.2 990.4
Coefficients:
Estimate Std. Error df t value Pr(>|t|)
period8 152.167 2.060 30.272 73.86 <0.0000000000000002 ***
period15 201.767 2.325 29.912 86.80 <0.0000000000000002 ***
period22 245.033 2.769 29.641 88.50 <0.0000000000000002 ***
period29 289.500 3.442 30.173 84.11 <0.0000000000000002 ***
period36 324.800 3.528 30.685 92.06 <0.0000000000000002 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Covariance estimate:
8 15 22 29 36
8 127.3458 133.0461 146.7499 168.9382 160.3436
15 133.0461 162.0991 178.7955 205.8290 195.3576
22 146.7499 178.7955 229.9818 264.7546 251.2854
29 168.9382 205.8290 264.7546 355.4306 337.3482
36 160.3436 195.3576 251.2854 337.3482 373.3904cov2cor(VarCorr(fit_ar1h)) 8 15 22 29 36
8 1.0000000 0.9260180 0.8575093 0.7940691 0.7353223
15 0.9260180 1.0000000 0.9260180 0.8575093 0.7940691
22 0.8575093 0.9260180 1.0000000 0.9260180 0.8575093
29 0.7940691 0.8575093 0.9260180 1.0000000 0.9260180
36 0.7353223 0.7940691 0.8575093 0.9260180 1.0000000model comparison
data.frame(model=c("UN","CS","CSH","ANTE"),
AIC=round(c(AIC(fit_full),AIC(fit_cs),AIC(fit_csh),AIC(fit_ante1)),3),
BIC=round(c(BIC(fit_full),BIC(fit_cs),BIC(fit_csh),BIC(fit_ante1)),3)) model AIC BIC
1 UN 1008.781 1029.799
2 CS 1096.548 1099.351
3 CSH 1067.045 1075.452
4 ANTE 1006.587 1019.198data.frame(model=c("UN","AR1H","TOEP","ANTE"),
AIC=round(c(AIC(fit_full),AIC(fit_ar1h),AIC(fit_toep),AIC(fit_ante1)),3),
BIC=round(c(BIC(fit_full),BIC(fit_ar1h),BIC(fit_toep),BIC(fit_ante1)),3)) model AIC BIC
1 UN 1008.781 1029.799
2 AR1H 1002.381 1010.788
3 TOEP 999.736 1012.347
4 ANTE 1006.587 1019.198Let’s test now the mean trend
fit_ante1<-mmrm(formula=weight~period+ar1h(period|indiv),data = long_rats)
summary(fit_ante1)mmrm fit
Formula: weight ~ period + ar1h(period | indiv)
Data: long_rats (used 150 observations from 30 subjects with maximum 5
timepoints)
Covariance: heterogeneous auto-regressive order one (6 variance parameters)
Method: Satterthwaite
Vcov Method: Asymptotic
Inference: REML
Model selection criteria:
AIC BIC logLik deviance
1002.4 1010.8 -495.2 990.4
Coefficients:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 152.1667 2.0603 30.2700 73.86 <0.0000000000000002 ***
period15 49.6000 0.8823 72.4200 56.22 <0.0000000000000002 ***
period22 92.8667 1.4586 62.6000 63.67 <0.0000000000000002 ***
period29 137.3333 2.1977 49.9100 62.49 <0.0000000000000002 ***
period36 172.6333 2.4498 47.2800 70.47 <0.0000000000000002 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Covariance estimate:
8 15 22 29 36
8 127.3458 133.0461 146.7499 168.9382 160.3436
15 133.0461 162.0991 178.7955 205.8290 195.3576
22 146.7499 178.7955 229.9818 264.7546 251.2854
29 168.9382 205.8290 264.7546 355.4306 337.3482
36 160.3436 195.3576 251.2854 337.3482 373.3904emmeans(fit_ante1,specs = "period") period emmean SE df lower.CL upper.CL
8 152 2.06 30.3 148 156
15 202 2.32 29.9 197 207
22 245 2.77 29.6 239 251
29 290 3.44 30.2 282 297
36 325 3.53 30.7 318 332
Confidence level used: 0.95 emmeans(fit_ante1,specs = "period")%>%joint_tests() #this is the only way to do an ANOVA with mmrm model term df1 df2 F.ratio p.value
period 4 47.28 1460.461 <.0001Anova(fit_ante1)Analysis of Fixed Effect Table (Type II F tests)
Num Df Denom Df F Statistic Pr(>=F)
period 4 83.262 1460.5 < 0.00000000000000022 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1emmeans(fit_ante1,specs = "period")%>%pairs() contrast estimate SE df t.ratio p.value
period8 - period15 -49.6 0.882 72.4 -56.218 <.0001
period8 - period22 -92.9 1.460 62.6 -63.667 <.0001
period8 - period29 -137.3 2.200 49.9 -62.489 <.0001
period8 - period36 -172.6 2.450 47.3 -70.468 <.0001
period15 - period22 -43.3 1.070 60.8 -40.352 <.0001
period15 - period29 -87.7 1.880 48.3 -46.702 <.0001
period15 - period36 -123.0 2.200 46.5 -56.006 <.0001
period22 - period29 -44.5 1.370 49.2 -32.574 <.0001
period22 - period36 -79.8 1.830 49.5 -43.516 <.0001
period29 - period36 -35.3 1.340 71.6 -26.281 <.0001
P value adjustment: tukey method for comparing a family of 5 estimates emmeans(fit_ante1,specs = "period")%>%contrast(.,method="consec") contrast estimate SE df t.ratio p.value
period15 - period8 49.6 0.882 72.4 56.218 <.0001
period22 - period15 43.3 1.070 60.8 40.352 <.0001
period29 - period22 44.5 1.370 49.2 32.574 <.0001
period36 - period29 35.3 1.340 71.6 26.281 <.0001
P value adjustment: mvt method for 4 tests emmeans(fit_ante1,specs = "period")%>%contrast(.,method="poly") contrast estimate SE df t.ratio p.value
linear 433.00 6.46 49.8 67.068 <.0001
quadratic -27.40 3.59 99.5 -7.639 <.0001
cubic -2.83 2.44 60.1 -1.159 0.2508
quartic -17.90 5.09 95.7 -3.519 0.0007Saturated mean and random coefficiets
long_rats<-mutate(long_rats,daysc=scale(days,center = T,scale = T),daysc2=daysc^2,daysc3=daysc^3)
long_rats# A tibble: 150 × 7
indiv days weight period daysc[,1] daysc2[,1] daysc3[,1]
<fct> <dbl> <int> <fct> <dbl> <dbl> <dbl>
1 1 8 151 8 -1.41 1.99 -2.80
2 1 15 199 15 -0.705 0.497 -0.350
3 1 22 246 22 0 0 0
4 1 29 283 29 0.705 0.497 0.350
5 1 36 320 36 1.41 1.99 2.80
6 2 8 145 8 -1.41 1.99 -2.80
7 2 15 199 15 -0.705 0.497 -0.350
8 2 22 249 22 0 0 0
9 2 29 293 29 0.705 0.497 0.350
10 2 36 354 36 1.41 1.99 2.80
# ℹ 140 more rowsfit_poly0<-lmer(formula=weight~period-1+(1|indiv),data = long_rats)
fit_poly1<-lmer(formula=weight~period-1+(1+daysc|indiv),data = long_rats)
fit_icm1<-lmer(formula=weight~period-1+(1|indiv)+(daysc-1|indiv),data = long_rats)
fit_poly2<-lmer(formula=weight~period-1+(1+daysc+daysc2|indiv),data = long_rats)
fit_poly3<-lmer(formula=weight~period-1+(1+daysc+daysc2+daysc3|indiv),data = long_rats)
print(VarCorr(fit_poly0),comp="Variance") Groups Name Variance
indiv (Intercept) 199.643
Residual 53.783 summary(fit_cs)mmrm fit
Formula: weight ~ period - 1 + cs(period | indiv)
Data: long_rats (used 150 observations from 30 subjects with maximum 5
timepoints)
Covariance: compound symmetry (2 variance parameters)
Method: Satterthwaite
Vcov Method: Asymptotic
Inference: REML
Model selection criteria:
AIC BIC logLik deviance
1096.5 1099.4 -546.3 1092.5
Coefficients:
Estimate Std. Error df t value Pr(>|t|)
period8 152.167 2.906 41.640 52.35 <0.0000000000000002 ***
period15 201.767 2.906 41.640 69.42 <0.0000000000000002 ***
period22 245.033 2.906 41.640 84.31 <0.0000000000000002 ***
period29 289.500 2.906 41.640 99.61 <0.0000000000000002 ***
period36 324.800 2.906 41.640 111.75 <0.0000000000000002 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Covariance estimate:
8 15 22 29 36
8 253.4262 199.6432 199.6432 199.6432 199.6432
15 199.6432 253.4262 199.6432 199.6432 199.6432
22 199.6432 199.6432 253.4262 199.6432 199.6432
29 199.6432 199.6432 199.6432 253.4262 199.6432
36 199.6432 199.6432 199.6432 199.6432 253.4262199.643 +53.783 [1] 253.426data.frame(model=c("RCM0","RCM1","ICM1","RCM2","RCM3"),
AIC=round(c(AIC(fit_poly0),AIC(fit_poly1),AIC(fit_icm1),AIC(fit_poly2),AIC(fit_poly3)),3),
BIC=round(c(BIC(fit_poly0),BIC(fit_poly1),BIC(fit_icm1),BIC(fit_poly2),BIC(fit_poly3)),3)) model AIC BIC
1 RCM0 1106.548 1127.623
2 RCM1 1031.183 1058.279
3 ICM1 1041.051 1065.136
4 RCM2 1020.587 1056.715
5 RCM3 1015.622 1063.792data.frame(model=c("UN","CS","CSH","ANTE"),
AIC=round(c(AIC(fit_full),AIC(fit_cs),AIC(fit_csh),AIC(fit_ante1)),3),
BIC=round(c(BIC(fit_full),BIC(fit_cs),BIC(fit_csh),BIC(fit_ante1)),3)) model AIC BIC
1 UN 1008.781 1029.799
2 CS 1096.548 1099.351
3 CSH 1067.045 1075.452
4 ANTE 1002.381 1010.788Compare CS models!
Data input 2
Second dataset: weights from 27 rats over several weeks, subject to three diets
library(nparLD)Loading required package: MASS
Attaching package: 'MASS'The following object is masked from 'package:dplyr':
selectdata(rat)
glimpse(rat)Rows: 135
Columns: 4
$ resp <int> 46, 49, 49, 51, 52, 56, 57, 57, 60, 63, 52, 52, 54, 56, 57, 59…
$ time <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ group <fct> control, control, control, control, control, control, control,…
$ subject <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,… rat<-mutate(rat,period=as.factor(time),animal=as.factor(subject),resp=as.double(resp),trt=group)%>%dplyr::select(resp,time,animal,trt, period)
glimpse(rat)Rows: 135
Columns: 5
$ resp <dbl> 46, 49, 49, 51, 52, 56, 57, 57, 60, 63, 52, 52, 54, 56, 57, 59,…
$ time <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ animal <fct> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, …
$ trt <fct> control, control, control, control, control, control, control, …
$ period <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …ggplot(rat,aes(x=time,y=resp,color=animal))+geom_point()+geom_line()+geom_smooth(aes(color=NULL))`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
ggplot(rat,aes(x=time,y=resp,group=animal,color=trt))+geom_point()+geom_line()+geom_smooth(aes(group=trt,color=trt),lwd=2)`geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Compare models
First we run a models with different repeated measures structures
fit_cs<-mmrm(formula=resp~period*trt+cs(period|animal),data = rat)
fit_ar1<-mmrm(formula=resp~period*trt+ar1(period|animal),data = rat)
fit_csh<-mmrm(formula=resp~period*trt+csh(period|animal),data = rat)
fit_arh1<-mmrm(formula=resp~period*trt+ar1h(period|animal),data = rat)
fit_toep<-mmrm(formula=resp~period*trt+toeph(period|animal),data = rat)Warning in mmrm(formula = resp ~ period * trt + toeph(period | animal), :
Divergence with optimizer L-BFGS-B due to problems: L-BFGS-B needs finite
values of 'fn'fit_ante<-mmrm(formula=resp~period*trt+adh(period|animal),data = rat)
fit_un<-mmrm(formula=resp~period*trt+us(period|animal),data = rat)
data.frame(model=c("CS","AR1","CSH","ARH1","TOEP","ANTE","UN"),
AIC=round(c(AIC(fit_cs),AIC(fit_ar1),AIC(fit_csh),AIC(fit_arh1),AIC(fit_toep),AIC(fit_ante),AIC(fit_un)),3),
BIC=round(c(BIC(fit_cs),BIC(fit_ar1),BIC(fit_csh),BIC(fit_arh1),BIC(fit_toep),BIC(fit_ante),BIC(fit_un)),3)) model AIC BIC
1 CS 901.113 903.705
2 AR1 823.634 826.226
3 CSH 850.132 857.907
4 ARH1 781.151 788.926
5 TOEP 771.296 782.959
6 ANTE 783.104 794.766
7 UN 768.100 787.538Compare means
fit_un<-mmrm(formula=resp~period*trt+us(period|animal),data = rat)
Anova(fit_un,type="III")Analysis of Fixed Effect Table (Type III F tests)
Num Df Denom Df F Statistic Pr(>=F)
period 4 24 334.10 < 0.00000000000000022 ***
trt 2 24 7.73 0.0025559 **
period:trt 8 24 6.68 0.0001304 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1trt_at_time<-emmeans(fit_un,~trt|period)
time_at_trt<-emmeans(fit_un,~period|trt)
trt_at_timeperiod = 0:
trt emmean SE df lower.CL upper.CL
control 54.0 1.47 24 51.0 57.0
thiour 54.7 1.47 24 51.7 57.7
thyrox 55.6 1.76 24 51.9 59.2
period = 1:
trt emmean SE df lower.CL upper.CL
control 78.5 2.62 24 73.1 83.9
thiour 76.3 2.62 24 70.9 81.7
thyrox 75.9 3.13 24 69.4 82.3
period = 2:
trt emmean SE df lower.CL upper.CL
control 106.0 3.08 24 99.6 112.4
thiour 95.8 3.08 24 89.4 102.2
thyrox 104.9 3.68 24 97.3 112.5
period = 3:
trt emmean SE df lower.CL upper.CL
control 130.1 4.26 24 121.3 138.9
thiour 108.2 4.26 24 99.4 117.0
thyrox 134.1 5.09 24 123.6 144.7
period = 4:
trt emmean SE df lower.CL upper.CL
control 160.6 5.18 24 149.9 171.3
thiour 124.0 5.18 24 113.3 134.7
thyrox 164.3 6.19 24 151.5 177.1
Confidence level used: 0.95 trt_at_time%>%pairs()period = 0:
contrast estimate SE df t.ratio p.value
control - thiour -0.700 2.08 24 -0.337 0.9395
control - thyrox -1.571 2.29 24 -0.686 0.7735
thiour - thyrox -0.871 2.29 24 -0.381 0.9235
period = 1:
contrast estimate SE df t.ratio p.value
control - thiour 2.200 3.71 24 0.593 0.8250
control - thyrox 2.643 4.09 24 0.647 0.7959
thiour - thyrox 0.443 4.09 24 0.108 0.9935
period = 2:
contrast estimate SE df t.ratio p.value
control - thiour 10.200 4.35 24 2.343 0.0688
control - thyrox 1.143 4.80 24 0.238 0.9692
thiour - thyrox -9.057 4.80 24 -1.888 0.1639
period = 3:
contrast estimate SE df t.ratio p.value
control - thiour 21.900 6.03 24 3.634 0.0036
control - thyrox -4.043 6.64 24 -0.609 0.8167
thiour - thyrox -25.943 6.64 24 -3.907 0.0019
period = 4:
contrast estimate SE df t.ratio p.value
control - thiour 36.600 7.33 24 4.995 0.0001
control - thyrox -3.686 8.07 24 -0.457 0.8920
thiour - thyrox -40.286 8.07 24 -4.990 0.0001
P value adjustment: tukey method for comparing a family of 3 estimates trt_at_time%>%plot()Error in ggplot2::theme_light(base_size = base_size, base_family = base_family, : unused arguments (header_family = header_family, ink = ink, paper = paper, accent = accent)time_at_trttrt = control:
period emmean SE df lower.CL upper.CL
0 54.0 1.47 24 51.0 57.0
1 78.5 2.62 24 73.1 83.9
2 106.0 3.08 24 99.6 112.4
3 130.1 4.26 24 121.3 138.9
4 160.6 5.18 24 149.9 171.3
trt = thiour:
period emmean SE df lower.CL upper.CL
0 54.7 1.47 24 51.7 57.7
1 76.3 2.62 24 70.9 81.7
2 95.8 3.08 24 89.4 102.2
3 108.2 4.26 24 99.4 117.0
4 124.0 5.18 24 113.3 134.7
trt = thyrox:
period emmean SE df lower.CL upper.CL
0 55.6 1.76 24 51.9 59.2
1 75.9 3.13 24 69.4 82.3
2 104.9 3.68 24 97.3 112.5
3 134.1 5.09 24 123.6 144.7
4 164.3 6.19 24 151.5 177.1
Confidence level used: 0.95 time_at_trt%>%plot()Error in ggplot2::theme_light(base_size = base_size, base_family = base_family, : unused arguments (header_family = header_family, ink = ink, paper = paper, accent = accent)time_at_trt%>%contrast(.,method="poly")trt = control:
contrast estimate SE df t.ratio p.value
linear 264.80 13.00 24 20.307 <.0001
quadratic 8.60 5.62 24 1.529 0.1393
cubic 3.40 3.30 24 1.030 0.3132
quartic 16.20 4.05 24 4.001 0.0005
trt = thiour:
contrast estimate SE df t.ratio p.value
linear 170.50 13.00 24 13.075 <.0001
quadratic -18.70 5.62 24 -3.325 0.0028
cubic 5.50 3.30 24 1.667 0.1086
quartic 15.50 4.05 24 3.828 0.0008
trt = thyrox:
contrast estimate SE df t.ratio p.value
linear 275.71 15.60 24 17.690 <.0001
quadratic 20.00 6.72 24 2.975 0.0066
cubic -7.86 3.94 24 -1.992 0.0579
quartic 9.00 4.84 24 1.860 0.0752