1
$\begingroup$

Parametric bootstraps are considered the gold standard to compute confidence intervals for linear mixed effects models. I think they also allow for easy re-transformation of model predictions to the original scale if the transformation was needed to meet the model assumptions.

Unfortunately, I have several doubts about the way I calculate p-values & effect sizes with PBs and apparently errors in my calculation.

The example data is taken from a BACI design experiment. Bee colonies were exposed to fields that were either treated with a pesticide or not. Virus levels before and after exposure were regarded to see whether there is a seasonal effect and whether the pesticide treatment affects the change over time.

I would like to know

  • how to calculate p-values and effect sizes
  • how to calculate estimates at each time point for each group with confidence intervals
  • whether to calculate the seasonal change based on both groups (temp) or only of the control group (tempA)
  • if the retransformation to the original scale is reasonable (in this case, but more so in a case when the presentation of results in the transformed scale is not convenient)

I am particularly concerned because the p-values differ if the model predictions were retransformed or not and because they differ strongly from p-values returned by the PBmodcomp function.

Data

mydata = structure(list(Treatment = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L), .Label = c("A", "B"), class = "factor"), 
    Period = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Before", 
    "After"), class = "factor"), Field_pair = structure(c(1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 4L, 4L, 4L, 4L, 
    4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 
    6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 2L, 2L, 2L, 2L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 
    5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 
    7L, 7L, 7L, 7L), .Label = c("P01", "P02", "P03", "P04", "P05", 
    "P10", "P11", "P12"), class = "factor"), Field = structure(c(12L, 
    12L, 12L, 12L, 6L, 6L, 6L, 6L, 1L, 1L, 1L, 1L, 16L, 16L, 
    16L, 16L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 
    11L, 11L, 11L, 11L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 12L, 
    12L, 12L, 12L, 6L, 6L, 6L, 6L, 1L, 1L, 1L, 1L, 16L, 16L, 
    16L, 16L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 
    11L, 11L, 11L, 11L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L), .Label = c("VR02", 
    "VR03", "VR04", "VR05", "VR06", "VR07", "VR09", "VR12", "VR13", 
    "VR14", "VR16", "VR17", "VR18", "VR20", "VR21", "VR23"), class = "factor"), 
    Colony = structure(c(8L, 47L, 53L, 73L, 37L, 38L, 49L, 74L, 
    12L, 52L, 55L, 79L, 21L, 34L, 84L, 92L, 14L, 51L, 77L, 81L, 
    1L, 35L, 75L, 96L, 6L, 27L, 91L, 95L, 3L, 15L, 18L, 87L, 
    22L, 25L, 78L, 82L, 43L, 56L, 72L, 90L, 8L, 47L, 53L, 73L, 
    37L, 38L, 49L, 74L, 12L, 52L, 55L, 79L, 21L, 34L, 84L, 92L, 
    14L, 51L, 77L, 81L, 1L, 35L, 75L, 96L, 6L, 27L, 91L, 95L, 
    3L, 15L, 18L, 87L, 22L, 25L, 78L, 82L, 43L, 56L, 72L, 90L
    ), .Label = c("1", "2", "3", "4", "5", "6", "7", "8", "9", 
    "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", 
    "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", 
    "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", 
    "40", "41", "42", "43", "44", "46", "47", "48", "49", "50", 
    "51", "52", "53", "54", "55", "56", "57", "72", "73", "74", 
    "75", "77", "79", "80", "81", "82", "83", "85", "86", "87", 
    "88", "90", "91", "92", "93", "94", "95", "96", "97", "98", 
    "99", "100", "101", "103", "104", "105", "106", "107", "108", 
    "109", "110", "111", "112", "113", "114", "115", "116"), class = "factor"), 
    Virus = c(64324.3824299998, 1827.14944545, 114029.587035, 
    810744.912599999, 594182.854649999, 22255.8398685, 18276.450108, 
    461371.340849999, 653.155131299999, 371277.560219999, 943.060861049998, 
    422965.026224999, 0, 1871.75032695, 532.237185899999, 1529314.6701, 
    0, 0, 1752.3190776, 56246.6672249999, 786957.775799999, 24.2034116939999, 
    0, 379553.501564999, 651172.869899999, 319788.320354999, 
    374994.300344999, 0, 696269.316749998, 3842.11815854999, 
    769117.423199999, 841.965529649997, 60409.4161649999, 655137.392699998, 
    1340999.8371, 321968.807894999, 363001.618874999, 141632.57703, 
    42593.8418324999, 273205.177454999, 142.475038125, 0, 0, 
    936.618511499998, 0, 67991.5660199998, 607.563119099998, 
    0, 0, 0, 876655.104149998, 283612.049804999, 0, 89.2017629999998, 
    0, 263640.7662, 0, 838.496572199998, 0, 0, 0, 0, 0, 290599.52124, 
    831.487359265629, 0, 7854.71079749999, 0, 174339.89013, 0, 
    0, 0, 35.4081442574999, 232.221923009999, 91.6300332149997, 
    730958.891249999, 0, 332524.349849999, 0, 0)), .Names = c("Treatment", 
"Period", "Field_pair", "Field", "Colony", "Virus"), row.names = c(NA, 
-80L), class = "data.frame")

Code:

library(lme4)
library(pbkrtest)
# Function to calculate wald-test p-values
p.value <- function(x) { 
    require(pbkrtest)
    # extract coefficients
    coefs <- data.frame(coef(summary(x)))
    # get the Kenward-Roger-approximated degrees of freedom                  
    df.KR <- get_ddf_Lb(x, fixef(x))
    # get p-values from the t-distribution using the t-values and approximated
    # degrees of freedom
    coefs$p.KR <- 2 * (1 - pt(abs(coefs$t.value), df.KR))
    coefs   
}

m1 = lmer(log(Virus+1) ~ Treatment * Period + (1|Field_pair/Field/Colony), data = mydata) 
summary(m1)
p.value(m1) # seasonal effect but no treatment effect
# no significant difference between the treatment groups to start with (p = .28)
# sig. seasonal effect (p < .0001)
# no sig. treatment effect (p = 0.76)

# Testing for treatment effect using LRT
m2 = update(m1, .~ Treatment + Period + (1|Field_pair/Field/Colony))

anova(m1, m2, test = "LRT") # p = 0.75

### ------
Virus_estimates = data.frame()
for(treat in c("A","B")){
    for(pd in c("Before","After")){
        newdat = subset(mydata, Treatment == treat & Period == pd)
        funcx <- function(.){
            newdat$Treatment = treat
            pred = predict(.,newdata = newdat)
            pred = mean(10^pred-1)
        }
        library(lme4)
        lmer_Virus = lmer(log10(Virus+1) ~ Treatment * Period + (1|Field_pair/Field/Colony), data = mydata)
        out_Virus <- bootMer(lmer_Virus, funcx, seed = 1234, nsim = 100)
        assign(paste("t","Virus",treat,pd, sep = "_"),out_Virus$t)
        Q = quantile(out_Virus$t, c(0.5,0.025,0.975))
        D = as.data.frame(out_Virus$t0)
        names(D) = "t0"
        D$median = as.numeric(Q[1])
        D$CI_lwr = as.numeric(Q[2])
        D$CI_upr = as.numeric(Q[3])
        D$Treatment = treat
        D$Period = pd
        Virus_estimates = rbind(Virus_estimates, D)
    }
}
Virus_estimates_log = data.frame(log10(Virus_estimates[,1:4]), Virus_estimates[,5:6])

# Same calculation without retransforming predictions

Virus_estimates2 = data.frame()
for(treat in c("A","B")){
    for(pd in c("Before","After")){
        newdat = subset(mydata, Treatment == treat & Period == pd)
        funcx <- function(.){
            newdat$Treatment = treat
            pred = predict(.,newdata = newdat)
            pred = mean(pred)
        }
        library(lme4)
        lmer_Virus = lmer(log10(Virus+1) ~ Treatment * Period + (1|Field_pair/Field/Colony), data = mydata)
        assign(paste("t","Virus2",treat, pd, sep = "_"),out_Virus$t)
        out_Virus <- bootMer(lmer_Virus, funcx, seed = 1234, nsim = 100)
        Q = quantile(out_Virus$t, c(0.5,0.025,0.975))
        D = as.data.frame(out_Virus$t0)
        names(D) = "t0"
        D$median = as.numeric(Q[1])
        D$CI_lwr = as.numeric(Q[2])
        D$CI_upr = as.numeric(Q[3])
        D$Treatment = treat
        D$Period = pd
        Virus_estimates2 = rbind(Virus_estimates2, D)
    }
}
Virus_estimates2
Virus_estimates_log 

# Calculating estimates
start_diff = t_Virus_A_Before - t_Virus_B_Before
quantile(start_diff, c(0.5,0.025,0.975))

temp_A = t_Virus_A_After - t_Virus_A_Before
quantile(temp_A, c(0.5,0.025,0.975)) # sig
temp_B = t_Virus_B_After - t_Virus_B_Before
quantile(temp_B, c(0.5,0.025,0.975)) 
temp = c(temp_A, temp_B)
quantile(temp , c(0.5,0.025,0.975)) 

temp_treat = temp_B - temp_A 
quantile(temp_treat , c(0.5,0.025,0.975)) 

## Calculating effect sizes and p-values ----
t0_Virus_A_Before   = subset(Virus_estimates, Treatment == "A" & Period == "Before")$t0
t0_Virus_A_After    = subset(Virus_estimates, Treatment == "A" & Period == "After")$t0
t0_temp_A           = t0_Virus_A_After - t0_Virus_A_Before

t0_Virus_B_Before   = subset(Virus_estimates, Treatment == "B" & Period == "Before")$t0
t0_Virus_B_After    = subset(Virus_estimates, Treatment == "B" & Period == "After")$t0
t0_temp_B           = t0_Virus_B_After - t0_Virus_B_Before

t0_temp_treat = t0_temp_B - t0_temp_A 

p_PB = function(t0, t){
    p_PB = mean(abs(t - mean(t)) >= abs(t0))
    p_PB
} 

p_PB(t0_temp_treat, temp_treat) # p-value for treatment effect on seasonal change = 0.07

PBmodcomp(m1, m2, seed = 123, nsim = 100) # LRT-p = 0.7211 # PBtest = 0.7327

## Calculating estimates
start_diff = t_Virus2_A_Before - t_Virus2_B_Before
quantile(start_diff, c(0.5,0.025,0.975))

temp2_A = t_Virus2_A_After - t_Virus2_A_Before
quantile(temp2_A, c(0.5,0.025,0.975)) 
temp2_B = t_Virus2_B_After - t_Virus2_B_Before
quantile(temp2_B, c(0.5,0.025,0.975)) 
temp2 = c(temp2_A, temp2_B)
quantile(temp2 , c(0.5,0.025,0.975)) 

temp2_treat = temp2_B - temp2_A 
quantile(temp2_treat , c(0.5,0.025,0.975)) 

## Calculating effect sizes and p-values ----
t0_Virus2_A_Before   = subset(Virus_estimates2, Treatment == "A" & Period == "Before")$t0
t0_Virus2_A_After    = subset(Virus_estimates2, Treatment == "A" & Period == "After")$t0
t0_temp2_A           = t0_Virus2_A_After - t0_Virus2_A_Before

t0_Virus2_B_Before   = subset(Virus_estimates2, Treatment == "B" & Period == "Before")$t0
t0_Virus2_B_After    = subset(Virus_estimates2, Treatment == "B" & Period == "After")$t0
t0_temp2_B           = t0_Virus2_B_After - t0_Virus2_B_Before

t0_temp2_treat = t0_temp2_B - t0_temp2_A 

p_PB = function(t0, t){
    p_PB = mean(abs(t - mean(t)) >= abs(t0))
    p_PB
} 

p_PB(t0_temp_treat, temp_treat) # p-value for treatment effect on seasonal change = 1

PBmodcomp(m1, m2, seed = 123, nsim = 100) # LRT-p = 0.7211 # PBtest = 0.7327

PS: I hope the long code didn't scare you off, half of it is just copy paste with a slight change. I asked a similar question on LRTs here.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.