> library("dplyr")

> library("lme4")

> library("performance")

> library("marginaleffects")

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

> d <- mutate(d, room = factor(room), calf = factor(calf), 
+     diet = factor(diet, levels = c("Control", "Enhanced")), time_factor = factor(time_we .... [TRUNCATED] 

> print(head(d, 3))
# A tibble: 3 × 6
  room  calf  diet    time_week weight_kg time_factor
  <fct> <fct> <fct>       <dbl>     <dbl> <fct>      
1 R1    C01   Control         0      47.9 0          
2 R1    C01   Control         2      58.6 2          
3 R1    C01   Control         4      65.8 4          

> print(count(distinct(d, calf, room, diet), room, diet))
# A tibble: 8 × 3
  room  diet         n
  <fct> <fct>    <int>
1 R1    Control      4
2 R1    Enhanced     4
3 R2    Control      4
4 R2    Enhanced     4
5 R3    Control      4
6 R3    Enhanced     4
7 R4    Control      4
8 R4    Enhanced     4

> plot(NA, xlim = c(0, 8), ylim = range(d$weight_kg), 
+     xlab = "Week", ylab = "Weight (kg)")

> for (id in levels(d$calf)) {
+     z <- filter(d, calf == id)
+     lines(z$time_week, z$weight_kg, col = ifelse(z$diet[1] == 
+         "Control",  .... [TRUNCATED] 

> legend("topleft", c("Control", "Enhanced"), col = c("#003d73", 
+     "#c54a20"), lty = 1, bty = "n")

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

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

REML criterion at convergence: 701.4

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.31755 -0.51062  0.00236  0.64098  1.95874 

Random effects:
 Groups   Name        Variance Std.Dev.
 calf     (Intercept) 29.094   5.394   
 Residual              2.541   1.594   
Number of obs: 160, groups:  calf, 32

Fixed effects:
                          Estimate Std. Error t value
(Intercept)               38.43207    2.18001  17.629
roomR2                     0.05863    2.72041   0.022
roomR3                     2.40960    2.72041   0.886
roomR4                     4.03348    2.72041   1.483
dietEnhanced               0.02500    1.98857   0.013
time_factor2              10.71919    0.56358  19.020
time_factor4              19.11194    0.56358  33.911
time_factor6              28.42313    0.56358  50.433
time_factor8              36.45800    0.56358  64.690
dietEnhanced:time_factor2  0.44562    0.79703   0.559
dietEnhanced:time_factor4  3.22037    0.79703   4.040
dietEnhanced:time_factor6  2.60100    0.79703   3.263
dietEnhanced:time_factor8  3.19725    0.79703   4.011

> print(check_convergence(calf_fit))
[1] TRUE
attr(,"gradient")
[1] 2.904132e-09

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

> calf_means <- avg_predictions(calf_fit, by = c("diet", 
+     "time_factor"), re.form = NA)

> print(calf_means)

     diet time_factor Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
 Control            0     40.1       1.41 28.5   <0.001 590.6  37.3   42.8
 Control            2     50.8       1.41 36.1   <0.001 946.1  48.0   53.5
 Control            4     59.2       1.41 42.1   <0.001   Inf  56.4   61.9
 Control            6     68.5       1.41 48.7   <0.001   Inf  65.7   71.2
 Control            8     76.5       1.41 54.4   <0.001   Inf  73.8   79.3
 Enhanced           0     40.1       1.41 28.5   <0.001 591.3  37.3   42.8
 Enhanced           2     51.2       1.41 36.4   <0.001 963.7  48.5   54.0
 Enhanced           4     62.4       1.41 44.4   <0.001   Inf  59.7   65.2
 Enhanced           6     71.1       1.41 50.6   <0.001   Inf  68.4   73.9
 Enhanced           8     79.7       1.41 56.7   <0.001   Inf  77.0   82.5

Type: response


> calf_gaps <- avg_comparisons(calf_fit, variables = "diet", 
+     by = "time_factor", re.form = NA)

> print(calf_gaps)

 time_factor Estimate Std. Error      z Pr(>|z|)   S  2.5 % 97.5 %
           0    0.025       1.99 0.0126    0.990 0.0 -3.872   3.92
           2    0.471       1.99 0.2367    0.813 0.3 -3.427   4.37
           4    3.245       1.99 1.6320    0.103 3.3 -0.652   7.14
           6    2.626       1.99 1.3206    0.187 2.4 -1.271   6.52
           8    3.222       1.99 1.6204    0.105 3.2 -0.675   7.12

Term: diet
Type: response
Comparison: Enhanced - Control


> print(check_model(calf_fit, check = c("linearity", 
+     "qq"), residual_type = "normal"))

> x <- readr::read_csv("data/cow-crossover.csv")

> x <- mutate(x, cow = factor(cow), sequence = factor(sequence, 
+     levels = c("AB", "BA")), period = factor(period, levels = c(1, 
+     2)), trea .... [TRUNCATED] 

> print(head(x, 4))
# A tibble: 4 × 5
  cow   sequence period treatment intake_kg_day
  <fct> <fct>    <fct>  <fct>             <dbl>
1 W01   AB       1      A                  21.6
2 W01   AB       2      B                  22.7
3 W02   BA       1      B                  22.5
4 W02   BA       2      A                  20.0

> print(count(distinct(x, cow, sequence), sequence))
# A tibble: 2 × 2
  sequence     n
  <fct>    <int>
1 AB          12
2 BA          12

> print(count(x, sequence, period, treatment))
# A tibble: 4 × 4
  sequence period treatment     n
  <fct>    <fct>  <fct>     <int>
1 AB       1      A            12
2 AB       2      B            12
3 BA       1      B            12
4 BA       2      A            12

> paired <- tidyr::pivot_wider(select(x, cow, treatment, 
+     intake_kg_day), names_from = treatment, values_from = intake_kg_day, 
+     names_pref .... [TRUNCATED] 

> plot(paired$intake_kg_day.A, paired$intake_kg_day.B, 
+     xlab = "Same cow under A (kg/day)", ylab = "Same cow under B (kg/day)")

> abline(a = 0, b = 1, lty = 2)

> crossover_fit <- lmer(intake_kg_day ~ treatment + 
+     period + sequence + (1 | cow), data = x, REML = TRUE)

> print(summary(crossover_fit))
Linear mixed model fit by REML ['lmerMod']
Formula: intake_kg_day ~ treatment + period + sequence + (1 | cow)
   Data: x

REML criterion at convergence: 180.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.70501 -0.37577  0.03845  0.41809  1.80745 

Random effects:
 Groups   Name        Variance Std.Dev.
 cow      (Intercept) 5.0542   2.2482  
 Residual             0.7033   0.8386  
Number of obs: 48, groups:  cow, 24

Fixed effects:
            Estimate Std. Error t value
(Intercept)  22.4265     0.6927  32.377
treatmentB    1.0401     0.2421   4.296
period2      -0.7238     0.2421  -2.990
sequenceBA   -2.2258     0.9492  -2.345

Correlation of Fixed Effects:
           (Intr) trtmnB perid2
treatmentB -0.175              
period2    -0.175  0.000       
sequenceBA -0.685  0.000  0.000

> print(check_convergence(crossover_fit))
[1] TRUE
attr(,"gradient")
[1] 6.234418e-08

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

> print(check_model(crossover_fit, check = c("linearity", 
+     "qq"), residual_type = "normal"))

> crossover_gap <- avg_comparisons(crossover_fit, variables = "treatment", 
+     re.form = NA)

> print(crossover_gap)

 Estimate Std. Error   z Pr(>|z|)    S 2.5 % 97.5 %
     1.04      0.242 4.3   <0.001 15.8 0.566   1.51

Term: treatment
Type: response
Comparison: B - A


> crossover_means <- avg_predictions(crossover_fit, 
+     by = "treatment", re.form = NA)

> print(crossover_means)

 treatment Estimate Std. Error    z Pr(>|z|)   S 2.5 % 97.5 %
         A       21       0.49 42.8   <0.001 Inf    20   21.9
         B       22       0.49 44.9   <0.001 Inf    21   23.0

Type: response