Skip to main content
3 of 3
added 188 characters in body
kjetil b halvorsen
  • 82.8k
  • 32
  • 201
  • 663

Multilevel-Regression Contrasts and Multiple Comparisons

My experiment looks roughly something like this: 5 Different biological replicates were tested with 4 Conditions (incl. Control) and 3 Time points (1h, 6h, 24h). My samples are nested within biological replicates (actually they are cells from primary neuron culture). Therefore, for each condition by time-point I have 15 (cells) per biological replicate.

Using Multilevel Regression I want to test whether there is a statistically significant difference between Conditions (4 levels) and Time points (3 levels) on a normalized score. A word on the outcome variable: The normalized score is the mean pixel intensity normalized to the respective control at the corresponding time point. Thus, the norm score 1.4 for condition x at time t would mean a 0.4 increase in mean pixel intensity compared to the control at time t.

I assume my question requires me to use either sum contrasts or effect coding followed by pairwise comparisons on the emms (lots of uncertainty here!).

My Questions:

  1. Is my approach reasonable or could I in fact directly test custom contrasts and perhaps save some multiple comparison correction? Which contrasts are best to use in this scenario?
  2. More general, the ICC in my null-model is below 5%. However, when I account for the predictor variables it increases. Which one is now correct and should be reported? Both?

I appreciate any additional comments to improve! This is my first Question. So, sorry in advance in case I have made an error in posing this problem!

Here is my example data frame with the initial analysis steps I have done. Please note that the example data in the Norm column is randomly generated and thus, significances will be different. In my case M3 was the "best" model and is thus of interest for me.

library(tidyverse)
library(lme4)

# Define levels for Experiment, Condition, and Time Point
experiment_levels <- c("Experiment1", "Experiment2", "Experiment3", 
                       "Experiment4", "Experiment5")
condition_levels <- c("Control", "NoSerum", "Rap", "EBSS")
time_levels <- c("1h", "6h", "24h")

# Number of observations per time point and condition
num_obs_per_condition <- 15

# Create data
df <- data.frame(
  Experiment = factor(rep(experiment_levels, 
          each = num_obs_per_condition * length(condition_levels) * 
                 length(time_point_levels))),
  Condition = factor(rep(rep(condition_levels, 
          each = num_obs_per_condition * length(time_point_levels)), 
  times = length(experiment_levels))),
  Time = factor(rep(rep(time_point_levels, 
          each = num_obs_per_condition), times = 
             length(experiment_levels) * length(condition_levels))), 
  Sample = rep(1:num_obs_per_condition, times = 
            length(experiment_levels) * length(condition_levels) * 
                        length(time_levels)),
  Norm = runif(num_obs_per_condition * length(experiment_levels) * 
    length(condition_levels) * length(time_point_levels), min = -2, 
           max = 2)
) %>% 
  mutate(Experiment = factor(Experiment, levels = experiment_levels),
         Condition = factor(Condition, levels = condition_levels),
         Time = factor(Time, levels = time_levels))

# Null Model vs Regression
reg <- lm(Norm~ 1, data = df)
M0 <- lmer(Norm ~ 1 + (1|Experiment), data = df, REML = FALSE)

# Log Likelihood Test to compare linear regression to random effects model
log1 <- logLik(reg)
log2 <- logLik(M0)
logRes <- 2*(log2-log1) 

# Model with Random Intercept
M1 <- lmer(Norm~ Condition * Time + (1|Experiment), REML = FALSE, 
           data = df)

# Model with RI+RS for Time
M2 <- lmer(Norm~ Condition * Time + (1 + Time|Experiment), 
           REML = FALSE, data = df)

# Model with RI + RS for Time and Condition
M3 <- lmer(Norm~Condition * Time + (1 + Time + Condition|Experiment), 
           REML = FALSE, data = df)

# Model with RI + RS for Time*Condition
M4 <- lmer(Norm~Condition * Time + (1 + Time * Condition|Experiment), 
           REML = FALSE, data = df)

# Test best Model
anova(M0, M1, M2, M3, M4, test = "Chisq")

# M3 with sum contrasts
contrasts(df$Condition) <- contr.sum(4)
contrasts(df$Time) <- contr.sum(3)
M3 <- lmer(Norm~Condition * Time + (1 + Time + Condition|Experiment), 
           REML = FALSE, data = df)
summary(M3)