• Statistical analysis of linear mixed-effects models is performed using the nlme package in R.
  • The reference book is “Pinheiro JC, Bates DM (2000). Mixed-Effects Models in S and S-PLUS. Springer, New York”.
  • The ggplot2 package in R is used to produce plots.
  • Please, install and load the packages.
# install.packages("nlme")  # Install the nlme package
library(nlme)               # Load the nlme package
library(tidyverse)            # Load the ggplot2 package

Section 1: A simple example of random effects

Objective

  • The dataset TumVol contains data from an experiment validating a new PDX model of non-small-cell lung cancer. Tumour cells were inoculated into six C57BL/6J mice.
  • Tumour volume (mm3) was measured after 2 weeks using a digital caliper.
  • Please, load and display the TumVol dataset as follows:
TumVol = read.csv("data/TumVol.csv") # Load the TumVol dataset
TumVol                               # Display the dataset
##    IDmouse TVolume
## 1        4     305
## 2        4     306
## 3        4     299
## 4        3     796
## 5        3     805
## 6        3     784
## 7        1     906
## 8        1     898
## 9        1     897
## 10       6     670
## 11       6     677
## 12       6     694
## 13       2     393
## 14       2     395
## 15       2     383
## 16       5     407
## 17       5     392
## 18       5     397

  • Tumour volume measurements were performed by the same operator.

  • They were repeated 3 times in the same day with 45 min intervals.

  • The quantities the researchers were interested in estimating from this experiment were:

    1. the average tumour volume for a “typical” PDX model (expected tumour volume)
    2. the variation in average tumour volume among PDXs (the between-mouse variability)
    3. the variation in the measured tumour volume for a single PDX (the within-mouse variability).

Let’s display the data and identify the parameters to estimate:

TumVol %>% 
  mutate(IDmouse = factor(IDmouse)) %>%
  ggplot(aes(x = TVolume, y = IDmouse, colour = IDmouse)) +
  geom_point()

  • We can see that there is considerable variability in the mean tumour volume for different mice.
  • Overall the between-mouse variability is much greater than the within-mouse variability.

Statistical analysis

  • Data like the TumVol example can be analyzed either with a fixed-effects model or with a random-effects model.
  • The distinction between the two regression models is according to whether we wish to make inferences about those particular mice that were used in the experiment or to make inferences about the population from which these mice were drawn.
  • Assume the simple model: yij = \(\beta\) + \(\epsilon\)ij, i=1,..,6; j=1,2,3. The errors \(\epsilon\)ij are assumed to be independent, identically and normally distributed.

Let’s fit the single-mean model with R:

fitModel1 <- lm(TVolume ~ 1, data=TumVol)
summary(fitModel1)
## 
## Call:
## lm(formula = TVolume ~ 1, data = TumVol)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -279.0 -185.8  -39.5  215.0  328.0 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   578.00      54.64   10.58 6.75e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 231.8 on 17 degrees of freedom
  • We can see that residual error is very large. The mean tumour volume is statistically different from zero.

Let’s ask ourself if the statistical model is a good model looking at the residuals behavior.

TumVol$resid <- resid(fitModel1)
TumVol %>% 
  mutate(IDmouse = factor(IDmouse)) %>%
  ggplot(aes(x = resid, y = IDmouse, colour = IDmouse)) +
  geom_point() +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red")

  • The statistical regression model assumes that errors are independent and identically distributed.
  • But errors of a single mouse are clearly correlated.
  • The previous assumption is strongly contradicted by real data.
  • Ignoring the classification factor (i.e. mice) when modeling grouped data leads to an inflated estimate of the within-mouse variability.

Fixed-effects models

To develop a better model, let’s introduce IDmouse as fixed-effects predictor:

TumVol$IDmouse <- factor(TumVol$IDmouse, ordered = FALSE, levels = c(1,2,3,4,5,6)) # IDmouse as categorical variable
fitModel2 <- lm(TVolume ~ IDmouse, data=TumVol)
summary(fitModel2)
## 
## Call:
## lm(formula = TVolume ~ IDmouse, data = TumVol)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.0000  -4.0833  -0.3333   4.1667  13.6667 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  900.333      4.724  190.59  < 2e-16 ***
## IDmouse2    -510.000      6.681  -76.34  < 2e-16 ***
## IDmouse3    -105.333      6.681  -15.77 2.19e-09 ***
## IDmouse4    -597.000      6.681  -89.36  < 2e-16 ***
## IDmouse5    -501.667      6.681  -75.09  < 2e-16 ***
## IDmouse6    -220.000      6.681  -32.93 3.89e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.182 on 12 degrees of freedom
## Multiple R-squared:  0.9991, Adjusted R-squared:  0.9988 
## F-statistic:  2727 on 5 and 12 DF,  p-value: < 2.2e-16

Residual error is strongly reduced respect to fitModel1.

Please, check the residuals behavior:

TumVol$resid <- resid(fitModel2)
TumVol %>% 
  mutate(IDmouse = factor(IDmouse)) %>%
  ggplot(aes(x = resid, y = IDmouse, colour = IDmouse)) +
  geom_point() +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red")

  • Test if variances are equal among mice groups using bartlett.test test
bartlett.test(resid ~ IDmouse, data = TumVol) # Bartlett’s homoscedasticity test
## 
##  Bartlett test of homogeneity of variances
## 
## data:  resid by IDmouse
## Bartlett's K-squared = 3.0898, df = 5, p-value = 0.6861
  • Histogram of residuals
# Overall
hist(TumVol$resid)

  • Check normality of residuals using Q-Q plot
qqnorm(TumVol$resid); qqline(TumVol$resid, col="red") # Q-Q plot

  • Check normality of residuals using Shapiro-Wilk normality test
shapiro.test(TumVol$resid) # Shapiro-Wilk normality test
## 
##  Shapiro-Wilk normality test
## 
## data:  TumVol$resid
## W = 0.97982, p-value = 0.9481
  • fitModel2 is definitively better than fitModel1, as confirmed by the ANOVA test:
anova(fitModel1, fitModel2) # ANOVA test
## Analysis of Variance Table
## 
## Model 1: TVolume ~ 1
## Model 2: TVolume ~ IDmouse
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     17 913646                                  
## 2     12    803  5    912843 2727.2 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • Even though the fixed-effects model accounts for the mouse effects, it does not provide a useful representation of the mice data for the following reasons:
    1. The number of parameters increases linearly with the number of mice.
    2. It only models the specific sample of mice used in the PDX experiment, while the main interest is in the population of mice from which the sample was drawn.
    3. This model does not provide an estimate of the between-mouse variability, which is one of the central quantities of interest in the PDX experiment.

Random-effects models

  • To solve the previous limits, we use a random-effects model with mice as random effects.
fitModel3 <- lme(TVolume ~ 1, random = ~1|IDmouse, data=TumVol)
summary(fitModel3)
## Linear mixed-effects model fit by REML
##   Data: TumVol 
##       AIC      BIC    logLik
##   168.155 170.6547 -81.07752
## 
## Random effects:
##  Formula: ~1 | IDmouse
##         (Intercept) Residual
## StdDev:    246.6452 8.181959
## 
## Fixed effects:  TVolume ~ 1 
##             Value Std.Error DF  t-value p-value
## (Intercept)   578   100.711 12 5.739197   1e-04
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.33469635 -0.50715104 -0.03989601  0.50084009  1.67492787 
## 
## Number of Observations: 18
## Number of Groups: 6
  • This model answers our 3 questions:
    1. the expected tumour volume = 578
    2. the between-mouse variability = 246.6452
    3. the within-mouse variability = 8.181959
  • What about the tumour volume of each mouse?
    • In fitModel2 mean tumour volume of each mouse is a parameter.
    • In fitModel3 mean tumour volume of each mouse is not a parameter. Their values are predicted (guessed) by the fitted model.

Please, refer to the following R commands:

rndEffects <- ranef(fitModel3)
rndEffects <- rndEffects$'(Intercept)'
rndEffects
## [1]  322.2151 -187.5979  216.9204 -274.5660 -179.2676  102.2958
  • Could fitModel3 reasonably generate our data?

Please, check the distribution of residuals and random effects.

TumVol$resid <- resid(fitModel3)
qqnorm(TumVol$resid); qqline(TumVol$resid, col="red") # Q-Q plot 

shapiro.test(TumVol$resid) # Shapiro-Wilk normality test
## 
##  Shapiro-Wilk normality test
## 
## data:  TumVol$resid
## W = 0.97943, p-value = 0.944
bartlett.test(resid ~ IDmouse, data = TumVol) # Bartlett’s homoscedasticity test
## 
##  Bartlett test of homogeneity of variances
## 
## data:  resid by IDmouse
## Bartlett's K-squared = 3.0898, df = 5, p-value = 0.6861
qqnorm(rndEffects); qqline(rndEffects, col="red") # Q-Q plot 

shapiro.test(rndEffects) # Shapiro-Wilk normality test
## 
##  Shapiro-Wilk normality test
## 
## data:  rndEffects
## W = 0.89655, p-value = 0.354
summary(rndEffects) # Random effects mean around zero
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -274.57 -185.52  -38.49    0.00  188.26  322.22
  • In conclusion, fitModel3 is a reasonable model and answer our questions. It could be used as “reference” regression model of the PDX experiment.

Section 2: Tumour growth curve analysis

Objective

  • The dataset TGCdSet contains data from an experiment performed in March 2024.
  • A preclinical murine model of non-small-cell lung cancer was used.

Please, load and display the first lines of the TGCdSet dataset.

TGCdSet = read.csv("data/TGCdSet.csv") # Load the TGCdSet dataset
head(TGCdSet)               # Display the first lines of the dataset
##   IDmouse ExpGroup       Date time  TVolume
## 1   ID-28  Vehicle 2024-03-18    0 663.6319
## 2   ID-16  Vehicle 2024-03-18    0 542.0157
## 3   ID-22  Vehicle 2024-03-18    0 391.3890
## 4    ID-5  Vehicle 2024-03-18    0 234.8196
## 5   ID-12  Vehicle 2024-03-18    0 394.5716
## 6   ID-15  Vehicle 2024-03-18    0 130.3431
# Data management
TGCdSet$IDmouse <- factor(TGCdSet$IDmouse, ordered = FALSE) # Categorical variable
TGCdSet$ExpGroup <- factor(TGCdSet$ExpGroup, ordered = FALSE, levels = c("Vehicle", "Carboplatin","Gefitinib","Carboplatin + Gefitinib")) # Categorical variable
  • On 18th March 2024, thirty-five mice were randomised to four experimental groups:
    1. Vehicle
    2. Carboplatin
    3. Gefitinib
    4. Carboplatin + Gefitinib
  • The primary outcome was tumour volume (mm3).
  • It was measured every 2-4 days using a digital caliper.
  • The primary objective was to statistically demonstrate (hypothesis testing) the synergistic effect of gefitinib combined with carboplatin.
  • As secondary objective researchers were interested to estimate the mean tumour volume and its 95% CI in the ‘vehicle’ group after two weeks (time = 14).

Statistical analysis

Exploratory analysis

  • We examine the data to …
    1. detect imbalance of baseline tumour volume between experimental groups
    2. to understand the shape of tumor growth curves.
  • Imbalance between experimental groups must be considered when interpreting results.
  • They should be avoided because it is impossible to disentangle treatment effect from the effect due to different baseline tumour volumes.
  • Randomisation techniques such as blocking should be applied to solve this methodological problem.
### Graphical checks
dBaseSet = subset(TGCdSet, time == 0)

ggplot(dBaseSet, aes(x=ExpGroup, y=TVolume)) + 
  geom_jitter(position=position_jitter(0.1)) +
  stat_summary(aes(x=ExpGroup, y=TVolume), fun = mean, colour="red", geom="point")

### Formal tools: hypothesis testing
kruskal.test(TVolume ~ ExpGroup, data = dBaseSet) # Kruskal-Wallis Rank Sum test
## 
##  Kruskal-Wallis rank sum test
## 
## data:  TVolume by ExpGroup
## Kruskal-Wallis chi-squared = 0.068427, df = 3, p-value = 0.9953
anovaTest = aov(TVolume ~ ExpGroup, data = dBaseSet) # ANOVA test
summary(anovaTest)
##             Df Sum Sq Mean Sq F value Pr(>F)
## ExpGroup     3   2097     699   0.028  0.994
## Residuals   30 750051   25002

  • No systematic difference between baseline values of experimental groups was detected.

  • Plots are primary tools to begin the modeling steps, specifically to ‘guess’ a reasonable regression model.

p <- ggplot(data=TGCdSet, aes(x=time, y=TVolume, group=IDmouse)) +
  geom_line( linetype = "solid",color="black") + 
  stat_summary( fun = mean, colour="red", geom="line", group=1)

p + 
  facet_grid(. ~ ExpGroup)

  • Tumour growth curves follows an exponential growth. To apply linear models, we work on natural log scale.
TGCdSet$TVlog = log(TGCdSet$TVolume)

p <- ggplot(data=TGCdSet, aes(x=time, y=TVlog, group=IDmouse)) +
  geom_line( linetype = "solid",color="black") + 
  stat_summary( fun = mean, colour="red", geom="line", group=1)
p + 
  facet_grid(. ~ ExpGroup)

  • Comment:
    • the plot suggests that growth trajectories on the natural log scale appear mostly straight, with slope’s dependence on experimental group.
    • Baseline tumor volumes are clearly different between mice. Moreover the ’Carboplatin + gefitinib group has a sample mean (red line) very different from the other three groups.

Model development

Based on the previous exploratory analysis a reasonable starting model is the following:

fitModel0 <- lme(TVlog ~ time + time:ExpGroup, 
                 random= ~1|IDmouse,  
                 method="REML", data = TGCdSet)
  • Based on residuals behavior, we could be satisfied by the starting model or we could be interested to improve it.
  • For instance, the exploratory analysis suggested that a common intercept could be used for all the experimental groups.
  • Moreover, the growth trajectories are linear on natural log scale but the slope seems to be different between mice.
  • Let’s formally check these two modification of the starting model using likelihood ratio test (LRT) and information criteria (i.e. AIC and BIC).
### Checking different slope by mouse
fitModel1 <- lme(TVlog ~ time + time:ExpGroup, 
                         random= ~1+time|IDmouse,  
                         method="REML", data = TGCdSet)
anova(fitModel1,fitModel0)
##           Model df      AIC      BIC     logLik   Test  L.Ratio p-value
## fitModel1     1  9 163.7722 198.4161  -72.88611                        
## fitModel0     2  7 495.2349 522.1802 -240.61746 1 vs 2 335.4627  <.0001
### Checking different intercept by experimental group
fitModel0ml <- lme(TVlog ~ time + time:ExpGroup, 
                           random= ~1+time|IDmouse,  
                           method="ML", data = TGCdSet)
fitModel1ml <- lme(TVlog ~ ExpGroup + time + time:ExpGroup, 
                           random= ~1+time|IDmouse,  
                           method="ML", data = TGCdSet)
anova(fitModel1ml,fitModel0ml)
##             Model df      AIC      BIC    logLik   Test   L.Ratio p-value
## fitModel1ml     1 12 138.9231 185.2866 -57.46153                         
## fitModel0ml     2  9 133.3955 168.1682 -57.69776 1 vs 2 0.4724599  0.9249

  • Lower values are better for both AIC and BIC.

  • AIC favors more complex models, while BIC includes a penalty for the number of parameters estimated so tends to favor more simple models with fewer parameters.

  • LRT requires that the compared models are nested. A nested model is simply one that contains a subset of the predictor variables in the overall regression model.

  • A likelihood ratio test uses the following null and alternative hypotheses:

    • HNull: The full model and the nested model fit the data equally well. Thus, you should use the nested model.
    • HAlternative: The full model fits the data significantly better than the nested model. Thus, you should use the full model.

  • If the p-value of the test is below a certain significance level (e.g. 0.05), then we can reject the null hypothesis and conclude that the full model offers a significantly better fit.

  • Important note: because the nested models fitModel0ml and fitModel1ml differ in the specification of their fixed-effects, a LRT can be defined for maximum likelihood estimator only.

  • Comment: the LRT, and the AIC and BIC criteria, all strongly favor the more general model with a different slope by mouse.

  • Let’s formally check if a quadratic term for the time predictor could improve the updated model.

### Checking a quadratic term for the time variable, overall and by experimental group
TGCdSet$time2 = (TGCdSet$time)^2
fitModel0ml <- lme(TVlog ~ time + time:ExpGroup, 
                           random= ~1+time|IDmouse,  
                           method="ML", data = TGCdSet)
fitModel2ml <- lme(TVlog ~ time + time2 + time:ExpGroup, 
                           random= ~1+time|IDmouse,  
                           method="ML", data = TGCdSet)
fitModel3ml <- lme(TVlog ~ time + time2 + (time+time2):ExpGroup, 
                           random= ~1+time|IDmouse,  
                           method="ML", data = TGCdSet)
anova(fitModel2ml,fitModel0ml)
##             Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## fitModel2ml     1 10 132.5431 171.1795 -56.27157                        
## fitModel0ml     2  9 133.3955 168.1682 -57.69776 1 vs 2 2.852382  0.0912
anova(fitModel3ml,fitModel0ml)
##             Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## fitModel3ml     1 13 122.9272 173.1544 -48.46358                        
## fitModel0ml     2  9 133.3955 168.1682 -57.69776 1 vs 2 18.46837   0.001

  • Comment: a quadratic term dependent on the experimental group seems to improve the model.

  • Now, let’s check the residuals of our updated model.

fitModel3 <- lme(TVlog ~ time + time2 + (time+time2):ExpGroup, 
                         random= ~1+time|IDmouse,  
                         method="REML", data = TGCdSet)
TGCdSet$resid <- resid(fitModel3) # Raw residuals
TGCdSet$fitted <- fitted(fitModel3) # Fitted values
  • Important note:
    • remember that residuals should be independent, normally distributed with mean zero and equal variance.
      • We strongly recommend to check residuals behavior by single and combined fixed and random effects.
ggplot(data=TGCdSet, aes(x= fitted, y=resid), group=IDmouse) +
  geom_point() + 
  geom_hline(yintercept=0, linetype="dashed", color = "blue") +
  labs(title = "Residuals by fitted value", x = "Fitted values", y = "Residuals")

ggplot(data=TGCdSet, aes(x= time, y=resid), group=IDmouse) +
  geom_point() + 
  geom_hline(yintercept=0, linetype="dashed", color = "blue") +
  labs(title = "Residuals by time", x = "Time", y = "Residuals") +
  stat_summary(aes(y=resid, group=1), fun = mean, colour="red", geom="line", group=1)

ggplot(data=TGCdSet, aes(x= time, y=resid), group=IDmouse) +
  geom_point() + 
  geom_hline(yintercept=0, linetype="dashed", color = "blue") +
  labs(title = "Residuals by time and experimental group", x = "Time", y = "Residuals") +
  stat_summary(aes(y=resid, group=1), fun = mean, colour="red", geom="line", group=1) + facet_grid(. ~ ExpGroup)

bartlett.test(resid ~ ExpGroup, data = TGCdSet) # Bartlett’s homoscedasticity test
## 
##  Bartlett test of homogeneity of variances
## 
## data:  resid by ExpGroup
## Bartlett's K-squared = 154.48, df = 3, p-value < 2.2e-16
ggplot(data=TGCdSet, aes(x= IDmouse, y=resid), group=IDmouse) +
  geom_point() + 
  geom_hline(yintercept=0, linetype="dashed", color = "blue") + 
  labs(title = "Residuals by mouse", x = "Mouse", y = "Residuals") +
  stat_summary(aes(y=resid), fun = mean, colour="red", geom="point", size=3)

ggplot(data=TGCdSet, aes(x= IDmouse, y=resid), group=IDmouse) +
  geom_point() + 
  geom_hline(yintercept=0, linetype="dashed", color = "blue") + 
  labs(title = "Residuals by mouse and time", x = "Mouse", y = "Residuals") +
  stat_summary(aes(y=resid), fun = mean, colour="red", geom="point", size=3) + facet_grid(. ~ ExpGroup)

  • Comment:
    • based on the residuals behaviour we detected a dependence of the variance on experimental groups (e.g. carbo + gefinitib has a larger variability).
  • We could solve the contradiction to the model’s assumptions using a different variance for each experimental group.
### Checking a quadratic term for the time variable, overall and by experimental group
fitModel3 <- lme(TVlog ~ time + time2 + (time+time2):ExpGroup, 
                         random= ~1+time|IDmouse,  
                         method="REML", data = TGCdSet)
fitModel4 <- lme(TVlog ~ time + time2 + (time+time2):ExpGroup, 
                         random= ~1+time|IDmouse, weights = varIdent(form = ~1|ExpGroup),
                         method="REML", data = TGCdSet)
anova(fitModel4,fitModel3)
##           Model df       AIC      BIC    logLik   Test  L.Ratio p-value
## fitModel4     1 16  97.96489 159.3686 -32.98245                        
## fitModel3     2 13 210.08213 259.9726 -92.04106 1 vs 2 118.1172  <.0001
  • Random effects (i.e. intercept and slope by mouse) are normally distributed around zero (statistics and plots not shown).
  • Based on residuals and random effects behaviour we consider the following model as “reference” (gold standard) model:
fitModel <- lme(TVlog ~ time + time2 + (time+time2):ExpGroup, 
                        random= ~1+time|IDmouse, weights = varIdent(form = ~1|ExpGroup), 
                        method="REML", data = TGCdSet)

Statistical inference

  • Let’s answer to the primary objective (hypothesis testing) and secondary objective (estimation) of the experiment. First of all, look at the output of the model fitting:
summary(fitModel)
## Linear mixed-effects model fit by REML
##   Data: TGCdSet 
##        AIC      BIC    logLik
##   97.96489 159.3686 -32.98245
## 
## Random effects:
##  Formula: ~1 + time | IDmouse
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev     Corr  
## (Intercept) 0.65180647 (Intr)
## time        0.03432278 0.047 
## Residual    0.28689477       
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | ExpGroup 
##  Parameter estimates:
## Carboplatin + Gefitinib                 Vehicle               Gefitinib 
##               1.0000000               0.3393315               0.3273767 
##             Carboplatin 
##               0.4911883 
## Fixed effects:  TVlog ~ time + time2 + (time + time2):ExpGroup 
##                                           Value  Std.Error  DF  t-value p-value
## (Intercept)                            5.617703 0.11254595 309 49.91475  0.0000
## time                                   0.101572 0.01252585 309  8.10899  0.0000
## time2                                 -0.000851 0.00024468 309 -3.47715  0.0006
## time:ExpGroupCarboplatin              -0.085980 0.01789576 309 -4.80448  0.0000
## time:ExpGroupGefitinib                -0.011562 0.01865658 309 -0.61973  0.5359
## time:ExpGroupCarboplatin + Gefitinib  -0.115881 0.01990542 309 -5.82155  0.0000
## time2:ExpGroupCarboplatin              0.001916 0.00031812 309  6.02162  0.0000
## time2:ExpGroupGefitinib                0.000475 0.00032099 309  1.47953  0.1400
## time2:ExpGroupCarboplatin + Gefitinib  0.001027 0.00040815 309  2.51540  0.0124
##  Correlation: 
##                                       (Intr) time   time2  tm:EGC tm:EGG t:EG+G
## time                                   0.002                                   
## time2                                  0.032 -0.363                            
## time:ExpGroupCarboplatin              -0.009 -0.700  0.254                     
## time:ExpGroupGefitinib                -0.005 -0.671  0.244  0.470              
## time:ExpGroupCarboplatin + Gefitinib  -0.059 -0.629  0.227  0.441  0.423       
## time2:ExpGroupCarboplatin              0.010  0.280 -0.768 -0.388 -0.188 -0.176
## time2:ExpGroupGefitinib               -0.005  0.277 -0.762 -0.194 -0.336 -0.174
## time2:ExpGroupCarboplatin + Gefitinib  0.069  0.218 -0.597 -0.153 -0.147 -0.539
##                                       t2:EGC t2:EGG
## time                                               
## time2                                              
## time:ExpGroupCarboplatin                           
## time:ExpGroupGefitinib                             
## time:ExpGroupCarboplatin + Gefitinib               
## time2:ExpGroupCarboplatin                          
## time2:ExpGroupGefitinib                0.586       
## time2:ExpGroupCarboplatin + Gefitinib  0.461  0.456
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.06254066 -0.49759350  0.02523054  0.52343675  3.09817920 
## 
## Number of Observations: 352
## Number of Groups: 35
  • Let’s prepare the matrix to specify the contrasts:
mTest = matrix(0, nrow = 2, ncol = 9)
# Synergy on the linear term
mTest[1,4] <- -1
mTest[1,5] <- -1
mTest[1,6] <- 1
# Synergy on the quadratic term
mTest[2,7] <- -1
mTest[2,8] <- -1
mTest[2,9] <- 1
mTest
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,]    0    0    0   -1   -1    1    0    0    0
## [2,]    0    0    0    0    0    0   -1   -1    1
  • Let’s perform the statistical test:
anova(fitModel, L=mTest)
## F-test for linear combination(s)
##   time:ExpGroupCarboplatin time:ExpGroupGefitinib
## 1                       -1                     -1
## 2                        0                      0
##   time:ExpGroupCarboplatin + Gefitinib time2:ExpGroupCarboplatin
## 1                                    1                         0
## 2                                    0                        -1
##   time2:ExpGroupGefitinib time2:ExpGroupCarboplatin + Gefitinib
## 1                       0                                     0
## 2                      -1                                     1
##   numDF denDF F-value p-value
## 1     2   309 6.07838  0.0026
# If ANOVA test with a F statistic (the gold standard test) could not be performed because the "between-" and "within-" variability is not clearly defined, you have other two options:
# 1. LRT with maximum likelihood estimator.
# 2. Wald test with a chi-square statistic:
#    library(multcomp)
#    glhtOutput <- glht(fitModel, linfct=mTest)
#    summary(glhtOutput, test = Chisqtest()) # Wald test

  • Comment:

  • The interaction effect of ‘carboplatin + gefitinib’ group has been statistically detected.

  • The interaction effect is related to the quadratic term of the time, not to the linear term of the time. Try to prove this!!!

  • Now, we estimate the mean tumour volume and its 95% CI in the ‘vehicle’ group after 14 days.

# Step n.1) parameters estimates are extracted from the model
fixedEffects <- as.numeric(fixef(fitModel))

# Step n.2) the point "vehicle group after 14 days" is identified by a design matrix
designMatrix <- matrix(0,nrow = 1, ncol = 9)
designMatrix[1,1] = 1
designMatrix[1,2] = 14
designMatrix[1,3] = 14^2

# Step n.3) point estimates and variance are calculated on log scale
pointEst = designMatrix %*% fixedEffects # point estimate on log scale
mVarCov <- as.matrix(vcov(fitModel)) # variance-covariance matrix on log scale
VarEst = designMatrix %*% mVarCov %*% t(designMatrix) # variance estimate on log scale

# Step n.4) 95% CI are calculated on log scale
low95CI = pointEst - sqrt(VarEst) * qnorm(0.975) # lower 95% CI in log scale
upp95CI = pointEst + sqrt(VarEst) * qnorm(0.975) # upper 95% CI in log scale

# Step n.5) point estimates and 95% CI are calculated on natural scale
pointEst = exp(pointEst)
low95CI = exp(low95CI)
upp95CI = exp(upp95CI)

c(pointEst,low95CI,upp95CI) # print
## [1]  965.7991  652.4579 1429.6215

Sensitivity analysis

  • Very important consideration:
    • “All models are wrong” 1
    • “Truth is much too complicated to allow anything but approximations” 2
    • The “reference” model is certainly a wrong model.
    • We must confirm its results assessing alternative models.

For instance let’s consider a model in which residuals of a mouse are not independent but correlated in consecutive time points:

fitSens <- lme(TVlog ~ time + time2 + (time+time2):ExpGroup, 
                        random= ~1+time|IDmouse, weights = varIdent(form = ~1|ExpGroup), correlation = corCAR1(form = ~time|IDmouse),
                        method="REML", data = TGCdSet)
summary(fitSens)
## Linear mixed-effects model fit by REML
##   Data: TGCdSet 
##         AIC       BIC   logLik
##   -71.26109 -6.019674 52.63055
## 
## Random effects:
##  Formula: ~1 + time | IDmouse
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev     Corr  
## (Intercept) 0.49626937 (Intr)
## time        0.02356297 0.546 
## Residual    0.71454857       
## 
## Correlation Structure: Continuous AR(1)
##  Formula: ~time | IDmouse 
##  Parameter estimate(s):
##       Phi 
## 0.9753731 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | ExpGroup 
##  Parameter estimates:
## Carboplatin + Gefitinib                 Vehicle               Gefitinib 
##               1.0000000               0.4194924               0.4935014 
##             Carboplatin 
##               0.5260473 
## Fixed effects:  TVlog ~ time + time2 + (time + time2):ExpGroup 
##                                           Value  Std.Error  DF  t-value p-value
## (Intercept)                            5.596914 0.10955396 309 51.08820  0.0000
## time                                   0.105882 0.01185109 309  8.93440  0.0000
## time2                                 -0.001101 0.00040224 309 -2.73711  0.0066
## time:ExpGroupCarboplatin              -0.078788 0.01734891 309 -4.54140  0.0000
## time:ExpGroupGefitinib                -0.018851 0.01854577 309 -1.01645  0.3102
## time:ExpGroupCarboplatin + Gefitinib  -0.122957 0.02323795 309 -5.29120  0.0000
## time2:ExpGroupCarboplatin              0.001822 0.00052671 309  3.45859  0.0006
## time2:ExpGroupGefitinib                0.000895 0.00060566 309  1.47694  0.1407
## time2:ExpGroupCarboplatin + Gefitinib  0.001365 0.00068115 309  2.00369  0.0460
##  Correlation: 
##                                       (Intr) time   time2  tm:EGC tm:EGG t:EG+G
## time                                   0.116                                   
## time2                                  0.017 -0.654                            
## time:ExpGroupCarboplatin              -0.030 -0.677  0.448                     
## time:ExpGroupGefitinib                -0.028 -0.634  0.419  0.432              
## time:ExpGroupCarboplatin + Gefitinib  -0.096 -0.514  0.333  0.347  0.324       
## time2:ExpGroupCarboplatin              0.004  0.501 -0.763 -0.679 -0.320 -0.256
## time2:ExpGroupGefitinib                0.008  0.437 -0.664 -0.298 -0.673 -0.223
## time2:ExpGroupCarboplatin + Gefitinib  0.014  0.389 -0.590 -0.265 -0.248 -0.769
##                                       t2:EGC t2:EGG
## time                                               
## time2                                              
## time:ExpGroupCarboplatin                           
## time:ExpGroupGefitinib                             
## time:ExpGroupCarboplatin + Gefitinib               
## time2:ExpGroupCarboplatin                          
## time2:ExpGroupGefitinib                0.507       
## time2:ExpGroupCarboplatin + Gefitinib  0.451  0.392
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.34694428 -0.47685179  0.03337708  0.43165290  1.74719681 
## 
## Number of Observations: 352
## Number of Groups: 35

The LRT shows that this more complex model could be used instead of the “reference” model.

anova(fitSens,fitModel)
##          Model df       AIC       BIC    logLik   Test L.Ratio p-value
## fitSens      1 17 -71.26109  -6.01967  52.63055                       
## fitModel     2 16  97.96489 159.36858 -32.98245 1 vs 2 171.226  <.0001

Let’s perform again the primary statistical test on this alternative model:

anova(fitSens, L=mTest)
## F-test for linear combination(s)
##   time:ExpGroupCarboplatin time:ExpGroupGefitinib
## 1                       -1                     -1
## 2                        0                      0
##   time:ExpGroupCarboplatin + Gefitinib time2:ExpGroupCarboplatin
## 1                                    1                         0
## 2                                    0                        -1
##   time2:ExpGroupGefitinib time2:ExpGroupCarboplatin + Gefitinib
## 1                       0                                     0
## 2                      -1                                     1
##   numDF denDF  F-value p-value
## 1     2   309 5.424113  0.0048

  • This alternative model confirms the results of the “reference” model.

  • We prefer the “reference” model instead of this alternative model because it is simpler.

  • Note: the contrasts matrix of this model is the same of the “reference” model. It is not enough!!!!!

  • Ideally, every statistical model that could reasonably generate our data should be assessed. Broader is the sensitivity analysis we perform, larger is our confidence on the results of the “reference” statistical model.

  • References:

    • 1 George E.P.Box, Journal of the American Statistical Association, Vol.71,No.356.(Dec.,1976),pp.791-799
    • 2 John von Neumann, 1947

Section 3: Practicals

(i) A factorial design: using genotype to estimate treatment effect

  • A factorial design was used to to detect a treatment effect and its dependency on the genotype (WT vs KO). Mice were randomised at time 0 and last observations were collected after 15 days. Data are reported in the file Exercise 1.csv. The primary outcome is tumour volume (mm3).
  • Please, answer the following questions:
    1. define on which scale you prefer to analyse tumour volume.
    2. develop a “reference” linear mixed-effects model.
    3. detect the treatment effect and its dependency on the genotype (WT vs KO).
    4. estimate population mean of the control-wild type (Ctr-WT) group and its 95% CI 12 days after randomisation.
    5. develop an alternative regression model to use in a sensitivity analysis.

(ii) A meta-analysis of four preclinical experiments

  • The experiment reported at point (i) was performed in other three laboratories. Raw data were pooled together to better estimate efficacy of the active treatment. Data are reported in file Exercise 2.csv.
  • Please, answer the following questions:
    1. is treatment effect different between laboratories?
    2. could you detect and estimate a laboratory effect?

(iii) A novel immunotherapy compound at CRUK CI

  • A preclinical in vivo antitumor activity experiment was performed at Cancer Research UK Cambridge Institute. The activity of a novel immunotherapy compound was tested against a vehicle group. Tumour volume (mm3) was measured at 30, 60 and 120 days after randomisation. Data are reported in the file Exercise 3.csv.
  • Please, answer the following questions:
    1. develop a “reference” linear mixed-effects model.
    2. test and estimate the treatment effect at 60 days after randomisation.
    3. develop an alternative regression model to use in a sensitivity analysis.