1
$\begingroup$

Please don’t be shocked by the amount of code in this post, this problem is (hopefully) not as complex as it seems. This dataset (https://www.dropbox.com/sh/dfn6n9ipqocd1mh/AAD74pdCRw8RXVBjpwAl_Ajja?dl=0) contains 40 participants of the study I ran and I put it up, along with the code, to help people recreating my problem.

> head(dataset) #Ignore the column X
    X Subject      Dilemma Inevitable Personal SelfBeneficial Resp2
1 721      21    Bike Week          1        0              1     0
2 722      21 Shark Attack          1        0              0     0
3 723      21        Ebola          1        1              0     1
4 725      21         Nuke          0        1              0     1
5 727      21      Rowboat          0        1              1     1
6 729      21  Cinderblock          0        0              0     0


Description of study:

In my study I was interested in how different contextual factors (Inevitable, Personal & SelfBeneficial) in a moral dilemma influence the judgement of whether killing someone to save multiple others is morally acceptable (Resp2 = 1) or unacceptable (Resp2 = 0). Every participant had to answer to 24 dilemmas of which every dilemma had a personal (Personal = 0) and an impersonal (Personal = 1) variant, that means one sentence in the dilemma was different within the dilemma. This factor was varied between subjects, so when participant A saw the impersonal variant of the “Bike Week” dilemma then participant B saw the personal variant of the “Bike Week” dilemma. The Inevitable and SelfBeneficial factor were fixed to the dilemma; that means for example “Shark Week” was, for all participants, inevitable (Inevitable = 1) and not self beneficial (SelfBeneficial = 0). So you could say the factor Personal was varied within-dilemma whereas the factors Inevitable and SelfBeneficial were varied between-dilemmas. All factors were fully crossed with each other.

Analytical approach:

As I wanted to calculate the probabilities of stating that killing was morally acceptable I conducted a multilevel logistic regression, in other words a generalised mixed model with a logit function. For this I used the glmer function of the R package lme4 (version 1.1-7). Following recommendations by Barr et al. (2013) I wanted to include the maximal random effect structure to reduce the probability of Type I errors. The maximal random effect structure here would be a by-subject random intercept as well as by-subject random slopes for Inevitable, Personal and SelfBeneficial. Furthermore a by-dilemma random intercept and a by-dilemma random slope for Personal. We will omit the by-dilemma random slope for Personal for now though as it does not influence the problem.

The problem:

When running the model in glmer I get the following output:

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation)['glmerMod']
Family: binomial  ( logit )
Formula: Resp2 ~ SelfBeneficial + Personal + Inevitable + (1 + SelfBeneficial + Personal + Inevitable | Subject) + (1 | Dilemma)
Data: dataset
Control: glmerControl(optCtrl = list(maxfun = 5e+05), optimizer = "bobyqa”)

> Random effects:
 Groups  Name           Variance Std.Dev. Corr             
 Subject (Intercept)    2.79199  1.6709                    
         SelfBeneficial 0.73805  0.8591   -0.50            
         Personal       0.07412  0.2723    0.43  0.53      
         Inevitable     0.43871  0.6624   -0.01  0.10  0.37
 Dilemma (Intercept)    2.60176  1.6130                    
Number of obs: 960, groups:  Subject, 40; Dilemma, 24
Fixed effects:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)      0.5726     0.6526   0.878    0.380    
SelfBeneficial  -1.0209     0.6993  -1.460    0.144    
Personal         0.7037     0.1806   3.895  9.8e-05 ***
Inevitable      -0.6884     0.6924  -0.994    0.320    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

As you can see the factors SelfBeneficial and Inevitable are not significant, however they should be as their significance has been shown in previous studies (however using anova instead of glmms). Most importantly: When I calculate bootstrapped confidence intervals (cutting of the 2.5 and 97.5 percentile with resampling only on the subject level) all of the three factors are significant as expected.

Result from bootstrapping (100 iterations):

> [1] "(Intercept)"
 2.5% 97.5% 
 0.17  1.10 
[1] "SelfBeneficial"
 2.5% 97.5% 
-1.58 -0.76 
[1] "Personal"
 2.5% 97.5% 
 0.46  1.07 
[1] "Inevitable"
 2.5% 97.5% 
-1.13 -0.44 

The questions:

Why are the results as indicated by the p-values of the glmer output and the bootstrapped confidence intervals so fundamentally different?

Am I over-fitting the model by adding too many random effects? A model without the by-dilemma random intercept yields the three significant factors but I want to include the maximal random effect structure since this has an effect on variables on the subject-level (not included in this example).

Do you spot any fundamental mistake I did in my experimental design or the bootstrapping?

Note that this is not a problem of too little bootstrap iterations (I tried a lot more), bias introduced via bootstrapping (bias is at between .1 or .2 in bootstrap samples with more iterations) or the omission of the by-dilemma random slope for Personal.

The code:

Here I added the code I used so you can reproduce/recreate my problem. Sorry that this is probably not the most elegant way of doing this.

dataset <- read.csv("datasetproblem.csv")
library(lme4)

#The model
mm <- glmer(Resp2 ~ SelfBeneficial+Personal+Inevitable + (1 + SelfBeneficial+Personal+Inevitable | Subject) + (1|Dilemma), family = binomial("logit"), data= dataset, control=glmerControl(optCtrl=list(maxfun=500000), optimizer = "bobyqa"))
summary(mm)

#Function needed in the bootstrap-function to exclude models that did not converge
check <- function(model){
  em <- 0
  tryCatch(model, warning=function(err) {    
    mes <- substr(err$message,1,19) 
    em <<- ifelse(mes == "failure to converge","not converged","other warning")
  })
  if(em == 0){return(model)}
  if(em != 0){return(0)}  
}

sa <- unique(dataset$Subject)

#Actual bootstrapping function
bootcon <- function(Model,n.sim=500,n.iter=50000,verbose=F){
  frameall <- data.frame()  
  for(j in 1:n.sim){
    sam1 <- sample(sa,40,replace=T)  

    bdat <- data.frame()
    for(i in sam1){
      ret <- dataset[dataset$Subject == i,]
      bdat <- rbind(bdat,ret)
    }

    wi <-as.vector(as.character(summary(Model)$call)[2])    

    modl <- check(glmer(wi,family = binomial("logit"), data= bdat, control=glmerControl(optCtrl=list(maxfun=n.iter), optimizer = "bobyqa"), verbose=verbose))
    if(class(modl)[1] == "glmerMod"){  
      q <- as.data.frame(summary(modl)$coefficients)[1:2]
      coefficient <- rownames(q)
      rownames(q) <- NULL
      Number <- rep(j,nrow(q))

      frapart <- cbind(coefficient,q,Number)  
      frameall <- rbind(frameall,frapart)
      print(paste("Iteration no. ",j,". System time: ",Sys.time()))}

  }
  return(frameall)
}


#Running the bootstrapping function
output <- bootcon(mm,100)  



#Analyse bootstraps
for(i in unique(output$coefficient)){
      print(i)
      #print(round(mean(output$Estimate[output$coefficient == i]),2 ))
      print(round(quantile(output$Estimate[output$coefficient == i],probs = c(0.025,0.975)),2))
} 
$\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.