> library("dplyr")

> library("lme4")

> library("performance")

> library("marginaleffects", warn.conflicts = FALSE)

> library("ggplot2")

> d <- readr::read_csv("data/calf-growth.csv")

> d <- mutate(d, diet = factor(diet, levels = c("Control", 
+     "Enhanced")), room = factor(room), calf = factor(calf), week = time_week)

> d <- mutate(d, visit = factor(week, levels = c(0, 
+     2, 4, 6, 8)))

> d <- arrange(d, calf, week)

> print(head(select(d, "calf", "diet", "week", "weight_kg"), 
+     4))
# A tibble: 4 × 4
  calf  diet     week weight_kg
  <fct> <fct>   <dbl>     <dbl>
1 C01   Control     0      47.9
2 C01   Control     2      58.6
3 C01   Control     4      65.8
4 C01   Control     6      73.1

> print(count(distinct(d, calf, diet), diet))
# A tibble: 2 × 2
  diet         n
  <fct>    <int>
1 Control     16
2 Enhanced    16

> print(count(d, week))
# A tibble: 5 × 2
   week     n
  <dbl> <int>
1     0    32
2     2    32
3     4    32
4     6    32
5     8    32

> observed <- ggplot(d, aes(week, weight_kg, group = calf, 
+     colour = diet)) + geom_line(alpha = 0.4) + labs(x = "Week", 
+     y = "Weight (kg)" .... [TRUNCATED] 

> print(observed)

> linear <- lmer(weight_kg ~ diet * week + room + (1 | 
+     calf), data = d, REML = TRUE)

> curved <- lmer(weight_kg ~ diet * week + I(week^2) + 
+     room + (1 | calf), data = d, REML = TRUE)

> print(summary(curved))
Linear mixed model fit by REML ['lmerMod']
Formula: weight_kg ~ diet * week + I(week^2) + room + (1 | calf)
   Data: d

REML criterion at convergence: 725.7

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.27062 -0.60674  0.06788  0.62813  2.19214 

Random effects:
 Groups   Name        Variance Std.Dev.
 calf     (Intercept) 29.084   5.393   
 Residual              2.595   1.611   
Number of obs: 160, groups:  calf, 32

Fixed effects:
                  Estimate Std. Error t value
(Intercept)       38.39388    2.17105  17.684
dietEnhanced       0.20787    1.95706   0.106
week               5.38766    0.16500  32.652
I(week^2)         -0.10708    0.01903  -5.628
roomR2             0.05862    2.72041   0.022
roomR3             2.40960    2.72041   0.886
roomR4             4.03347    2.72041   1.483
dietEnhanced:week  0.42749    0.09005   4.747

Correlation of Fixed Effects:
            (Intr) dtEnhn week   I(w^2) roomR2 roomR3 roomR4
dietEnhancd -0.451                                          
week        -0.110  0.050                                   
I(week^2)    0.070  0.000 -0.923                            
roomR2      -0.627  0.000  0.000  0.000                     
roomR3      -0.627  0.000  0.000  0.000  0.500              
roomR4      -0.627  0.000  0.000  0.000  0.500  0.500       
dtEnhncd:wk  0.083 -0.184 -0.273  0.000  0.000  0.000  0.000

> linear_ml <- update(linear, REML = FALSE)

> curved_ml <- update(curved, REML = FALSE)

> print(AIC(linear_ml, curved_ml))
          df      AIC
linear_ml  9 774.4759
curved_ml 10 747.5685

> print(check_singularity(curved))
[1] FALSE

> print(check_normality(curved))
OK: residuals appear as normally distributed (p = 0.268).

> checks_linear <- check_model(linear, check = c("qq", 
+     "linearity", "homogeneity"))

> checks_curved <- check_model(curved, check = c("qq", 
+     "linearity", "homogeneity"))

> print(checks_linear)

> print(checks_curved)

> d <- mutate(d, error_kg = residuals(curved, type = "response"))

> time_check <- ggplot(d, aes(week, error_kg)) + geom_point(alpha = 0.35) + 
+     stat_summary(fun = mean, geom = "line", colour = "darkorange", 
+   .... [TRUNCATED] 

> print(time_check)

> grid <- mutate(tidyr::crossing(week = c(0, 4, 8), 
+     diet = d$diet, room = d$room), calf = d$calf[1])

> means <- avg_predictions(curved, newdata = grid, by = c("week", 
+     "diet"), re.form = NA)

> print(means)

 week     diet Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
    0 Control      40.0       1.39 28.7   <0.001 601.2  37.3   42.7
    0 Enhanced     40.2       1.39 28.9   <0.001 607.4  37.5   43.0
    4 Control      59.9       1.37 43.7   <0.001   Inf  57.2   62.5
    4 Enhanced     61.8       1.37 45.1   <0.001   Inf  59.1   64.5
    8 Control      76.3       1.39 54.8   <0.001   Inf  73.5   79.0
    8 Enhanced     79.9       1.39 57.4   <0.001   Inf  77.2   82.6

Type: response


> differences <- avg_comparisons(curved, newdata = grid, 
+     variables = list(diet = c("Control", "Enhanced")), by = "week", 
+     re.form = NA)

> print(differences)

 week Estimate Std. Error     z Pr(>|z|)   S  2.5 % 97.5 %
    0    0.208       1.96 0.106   0.9154 0.1 -3.628   4.04
    4    1.918       1.92 0.997   0.3188 1.6 -1.852   5.69
    8    3.628       1.96 1.854   0.0638 4.0 -0.208   7.46

Term: diet
Type: response
Comparison: Enhanced - Control


> mean_plot <- ggplot(means, aes(week, estimate, colour = diet, 
+     fill = diet)) + geom_ribbon(aes(ymin = conf.low, ymax = conf.high), 
+     alph .... [TRUNCATED] 

> print(mean_plot)

> d <- mutate(d, visit_index = match(week, c(0, 2, 4, 
+     6, 8)))

> library("nlme")

> ar1 <- nlme::lme(weight_kg ~ diet * week + I(week^2) + 
+     room, random = ~1 | calf, correlation = nlme::corAR1(form = ~visit_index | 
+     calf .... [TRUNCATED] 

> print(summary(ar1))
Linear mixed-effects model fit by REML
  Data: d 
       AIC      BIC    logLik
  737.8064 771.0691 -357.9032

Random effects:
 Formula: ~1 | calf
        (Intercept) Residual
StdDev:    5.324848  1.85135

Correlation Structure: AR(1)
 Formula: ~visit_index | calf 
 Parameter estimate(s):
      Phi 
0.4160279 
Fixed effects:  weight_kg ~ diet * week + I(week^2) + room 
                     Value Std.Error  DF   t-value p-value
(Intercept)       38.54645 2.1790112 125 17.689882  0.0000
dietEnhanced       0.19987 1.9768356  27  0.101104  0.9202
week               5.40378 0.1748254 125 30.909576  0.0000
I(week^2)         -0.10741 0.0194670 125 -5.517599  0.0000
roomR2            -0.10715 2.7224800  27 -0.039356  0.9689
roomR3             2.07494 2.7224800  27  0.762152  0.4526
roomR4             3.82653 2.7224800  27  1.405531  0.1713
dietEnhanced:week  0.41318 0.1123416 125  3.677918  0.0003
 Correlation: 
                  (Intr) dtEnhn week   I(w^2) roomR2 roomR3 roomR4
dietEnhanced      -0.454                                          
week              -0.116  0.073                                   
I(week^2)          0.056  0.000 -0.891                            
roomR2            -0.625  0.000  0.000  0.000                     
roomR3            -0.625  0.000  0.000  0.000  0.500              
roomR4            -0.625  0.000  0.000  0.000  0.500  0.500       
dietEnhanced:week  0.103 -0.227 -0.321  0.000  0.000  0.000  0.000

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.91591554 -0.52932423  0.08089803  0.56175374  1.99329556 

Number of Observations: 160
Number of Groups: 32 

> print(nlme::VarCorr(ar1))
calf = pdLogChol(1) 
            Variance  StdDev  
(Intercept) 28.354002 5.324848
Residual     3.427495 1.851350

> curved_ml <- update(curved, REML = FALSE)

> ar1_ml <- update(ar1, method = "ML")

> print(AIC(curved_ml, ar1_ml))
          df      AIC
curved_ml 10 747.5685
ar1_ml    11 740.5392

> week8 <- filter(grid, week == 8)

> extension_difference <- avg_comparisons(ar1, newdata = week8, 
+     variables = list(diet = c("Control", "Enhanced")), level = 0)

> print(extension_difference)

 Estimate Std. Error    z Pr(>|z|)   S  2.5 % 97.5 %
     3.51       1.98 1.77   0.0762 3.7 -0.369   7.38

Term: diet
Type: response
Comparison: Enhanced - Control