One-sided tests for Linear Mixed Models

Author

Lukas Graz

Published

November 29, 2024

  1. Simulate (and show) some Data
  2. fit linear mixed model
  3. anova and two-sided t-test
  4. one-sided t-test

Load Data

# read 20241202_HT_MC38-tumor volume_last day.xlsx file
library(readxl)
try(setwd("./posts/lmm-one-sided-t-tests"), silent=TRUE)

D <- read_excel("20241202_HT_MC38-tumor volume_last day.xlsx")
D$diameter <- D$V13 ^ (1/3) 
knitr::kable(D)
MiceID cage Sex Age group repetition V13 diameter
930 515365 M 13 KO 1 410.52384 7.432121
932 515365 M 13 KO 1 579.40350 8.336691
931 515365 M 13 WT 1 1329.90461 10.996982
938 515366 F 12 WT 1 830.65883 9.400282
939 515366 F 12 WT 1 996.99995 9.989990
933 515366 F 13 WT 1 372.07621 7.192457
934 515366 F 13 KO 1 421.96415 7.500528
935 516557 M 12 KO 1 776.17485 9.190092
936 516557 M 12 WT 1 1707.56800 11.952516
952 519341 M 9 WT 1 967.74602 9.891310
951 519341 M 9 KO 1 259.42441 6.377791
950 519341 M 9 WT 1 438.63516 7.598033
954 519342 F 9 WT 1 830.00991 9.397834
953 519342 F 9 WT 1 903.94049 9.668964
955 519342 F 9 WT 1 1604.01668 11.705850
888 518516 M 15 KO 2 849.40284 9.470464
887 518516 M 15 WT 2 518.76835 8.035098
894 520161 F 15 WT 2 306.81635 6.744651
890 520161 F 15 KO 2 341.42484 6.989268
895 521046 M 12 KO 2 293.99029 6.649326
896 521046 M 12 KO 2 686.81700 8.822947
897 521047 F 12 KO 2 503.68519 7.956457
900 521047 F 12 WT 2 518.97852 8.036183
898 521047 F 12 KO 2 35.89571 3.298736
899 521047 F 12 WT 2 241.72566 6.229324
901 523672 M 9 WT 2 348.53308 7.037439
903 523672 M 9 KO 2 328.50457 6.899969
902 523672 M 9 WT 2 312.02899 6.782633
931 524731 F 7 WT 2 1178.05030 10.561395
930 524731 F 7 WT 2 796.13798 9.268215
929 524731 F 7 WT 2 482.51191 7.843370
905 523673 F 11 WT 3 134.89640 5.128615
906 523673 F 11 KO 3 451.55379 7.671904
920 524714 M 9 KO 3 49.31999 3.667254
922 524714 M 9 WT 3 188.00473 5.728702
921 524714 M 9 KO 3 358.24065 7.102179
925 524715 M 9 WT 3 117.07240 4.891982
923 524715 M 9 KO 3 235.87591 6.178663
924 524715 M 9 WT 3 116.18825 4.879636
913 524737 M 9 WT 3 1127.79315 10.409019
907 524737 M 10 KO 3 304.16803 6.725189
908 524737 M 10 WT 3 505.34612 7.965193
914 524738 F 9 WT 3 1097.82540 10.315994
910 524738 F 10 WT 3 130.89710 5.077423
912 524738 F 10 WT 3 660.29680 8.707893
916 524739 F 9 KO 3 46.47149 3.595248
915 524739 F 9 WT 3 389.76231 7.304659
917 524739 F 9 WT 3 557.56208 8.230592

Since Volume is heavily right-skewed, we will use the cubic root of the volume as the response variable. The resulting variable approximately resembles the diameter.

fit linear mixed model

library(lmerTest)
D$group <- factor(D$group)

# if we only analyze the first repetition:
# fit <- lmer(diameter ~ group +        (1|cage), data=D)

# otherwise:
fit <- lmer(diameter ~ group + repetition + (1|cage), data=D)

anova and two-sided t-test

anova(fit)
Type III Analysis of Variance Table with Satterthwaite's method
           Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
group      17.386  17.386     1 42.403  6.5197 0.014352 * 
repetition 36.722  36.722     1 14.215 13.7702 0.002272 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(fit) # look at the p-value of groupWT
Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: diameter ~ group + repetition + (1 | cage)
   Data: D

REML criterion at convergence: 185

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.05849 -0.80721  0.05774  0.76071  1.86338 

Random effects:
 Groups   Name        Variance Std.Dev.
 cage     (Intercept) 0.2457   0.4957  
 Residual             2.6668   1.6330  
Number of obs: 48, groups:  cage, 17

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)   9.4218     0.7860 18.6525  11.987 3.32e-10 ***
groupWT       1.2757     0.4996 42.4025   2.553  0.01435 *  
repetition   -1.2116     0.3265 14.2151  -3.711  0.00227 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
           (Intr) gropWT
groupWT    -0.403       
repetition -0.854  0.010

So groupWT is 1.276 bigger than groupKO. The p-value is 0.014.

one-sided t-test

see this stackoverflow post

one_sided_ttest <- function(fit, which_term="group", greater_than_reference=TRUE){
  COEFS <- coef(summary(fit))
  stopifnot(any(grepl(which_term, names(fixef(fit)))))
  which_term <- grep(which_term, rownames(COEFS))
  t_val <- COEFS[which_term, "t value"]
  df <- COEFS[which_term, "df"]
  p_val <- pt(t_val, df, lower.tail=!greater_than_reference)
  cat("p-value for one-sided t-test for the alternative hypothethis:\n", 
    levels(D$group)[1], if(greater_than_reference) "<" else ">", levels(D$group)[2], ":", p_val, "\n")
  invisible(p_val)
}

# PVALUE: one-sided t-test for groupWT > groupKO
one_sided_ttest(fit, greater_than_reference=TRUE)
p-value for one-sided t-test for the alternative hypothethis:
 KO < WT : 0.007176008 
# PVALUE: one-sided t-test for groupWT < groupKO
one_sided_ttest(fit, greater_than_reference=FALSE)
p-value for one-sided t-test for the alternative hypothethis:
 KO > WT : 0.992824 

The one-sided p-value corresponds to half the two-sided p-value, if and only if the reference group is the one supposed to be smaller. The reference group is a term in R, and is (unless specified) the groupname that comes first in the alphabet. In this case, the reference group is groupKO. So the one-sided p-value for groupWT > groupKO is 0.007 which also equals 0.014/2.

Group Sequential Design

The consultancy case that has motivated this post, is doing a group sequential design. The idea is to stop the experiment early if the effect is very strong. The interim analysis is done after each repetition. The stopping rule is to stop if the one-sided p-value is smaller than 0,0221 (using the correction method described in Pocock (1977)).

# model with only the first repetition
fit1 <- lmer(diameter ~ group + (1|cage), 
             data=D, subset= repetition==1)
one_sided_ttest(fit1, greater_than_reference=TRUE)
p-value for one-sided t-test for the alternative hypothethis:
 KO < WT : 0.004332269 
# model with the first two repetitions
fit2 <- lmer(diameter ~ group + repetition + (1|cage), 
             data=D, subset= repetition %in% 1:2)
one_sided_ttest(fit2, greater_than_reference=TRUE)
p-value for one-sided t-test for the alternative hypothethis:
 KO < WT : 0.01616867 
# model with all three repetitions
fit3 <- lmer(diameter ~ group + repetition + (1|cage), 
             data=D, subset= repetition %in% 1:3)
one_sided_ttest(fit3, greater_than_reference=TRUE)
p-value for one-sided t-test for the alternative hypothethis:
 KO < WT : 0.007176008 

The first p-value is smaller than 0.0221, so we could have stopped the experiment after the first repetition.