How to calculate confidence intervals of lmer estimates using Parametric Bootstrap when variable has been transformed?

Question

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.