I have used the mixed_model function from GLMMadaptive to run a model with fixed effects and a random effect to account for clustering. I have extracted the ICC using icc() from performance. I want to estimate the ICC confidence intervals.

I know that other functions are available to do this that are compatible with lme4 models, but they don't seem to apply to GLMMadaptive models. (note that I have been advised to use GLMMadaptive rather than lme4 due to the differences in approximation methods.)

In the meantime, I have used simulate() to try to estimate ICC. This involved simulating new data from my model, re-fitting a new model to simulated data, and extracting parameters of interest. It looks like my code to update the model (newfit = update(m, newy ~ .) is not working and it produces warnings. When I explore the model output, the model output/estimates are the same as the original model and the ICC confidence interval output are just the ICC value.

In summary, my questions are:

- Is simulate() the right approach to estimate ICC confidence intervals from my model (or is there another specific function that will estimate the parameters)?
- If simulate() is the right approach, how can I fix my code to produce meaningful/correct estimates?

#Packages

library(magrittr)

library(dplyr)

#>

#> Attaching package: 'dplyr'

#> The following objects are masked from 'package:stats':

#>

#> filter, lag

#> The following objects are masked from 'package:base':

#>

#> intersect, setdiff, setequal, union

library(GLMMadaptive)

library(performance)

library(reprex)

#Create dataset

set.seed(432)

dummy_vec <- rbinom(400, 1, prob = 0.3) # Create dummy vector

head(dummy_vec) # Print dummy vector

#> [1] 0 0 1 0 1 1

dummy_mat <- matrix(dummy_vec, ncol = 4) # Convert vector to matrix

data <- as.data.frame(dummy_mat) # Convert matrix to data frame

cols<-c("2","3","4")

data %<>%

mutate_each_(funs(factor(.)),cols) # Convert columns to factors

#> Warning: `mutate_each_()`

was deprecated in dplyr 0.7.0.

#> Please use `across()`

instead.

#> This warning is displayed once every 8 hours.

#> Call `lifecycle::last_lifecycle_warnings()`

to see where this warning was generated.

#> Warning: `funs()`

was deprecated in dplyr 0.8.0.

#> Please use a list of either functions or lambdas:

#>

#> # Simple named list:

#> list(mean = mean, median = median)

#>

#> # Auto named with `tibble::lst()`

:

#> tibble::lst(mean, median)

#>

#> # Using lambdas

#> list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))

#> This warning is displayed once every 8 hours.

#> Call `lifecycle::last_lifecycle_warnings()`

to see where this warning was generated.

str(data)

#> 'data.frame': 100 obs. of 4 variables:

#> V1: int 0 0 1 0 1 1 1 0 0 1 ...
#> V2: Factor w/ 2 levels "0","1": 1 2 2 2 1 1 1 2 2 1 ...

#> V3: Factor w/ 2 levels "0","1": 1 2 1 2 2 1 2 2 1 2 ...
#> V4: Factor w/ 2 levels "0","1": 1 1 1 2 1 1 1 1 2 1 ...

names(data)[1] <- "distress" # Rename column names

names(data)[2] <- "mhhx"

names(data)[3] <- "gender"

names(data)[4] <- "exposure"

data$school<-sample(100, size = nrow(data), replace = TRUE) #Add new column

head(data) # Visualise data

#> distress mhhx gender exposure school

#> 1 0 0 0 0 4

#> 2 0 1 1 0 20

#> 3 1 1 0 0 21

#> 4 0 1 1 1 49

#> 5 1 0 1 0 94

#> 6 1 0 0 0 18

str(data)

#> 'data.frame': 100 obs. of 5 variables:

#> distress: int 0 0 1 0 1 1 1 0 0 1 ...
#> mhhx : Factor w/ 2 levels "0","1": 1 2 2 2 1 1 1 2 2 1 ...

#> gender : Factor w/ 2 levels "0","1": 1 2 1 2 2 1 2 2 1 2 ...
#> exposure: Factor w/ 2 levels "0","1": 1 1 1 2 1 1 1 1 2 1 ...

#> $ school : int 4 20 21 49 94 18 45 23 38 53 ...

#Run mixed model

m<-mixed_model(fixed=distress~mhhx+gender+exposure,random=~1|school,data=data, family="binomial")

summary(m)

#>

#> Call:

#> mixed_model(fixed = distress ~ mhhx + gender + exposure, random = ~1 |

#> school, data = data, family = "binomial")

#>

#> Data Descriptives:

#> Number of Observations: 100

#> Number of Groups: 67

#>

#> Model:

#> family: binomial

#> link: logit

#>

#> Fit statistics:

#> log.Lik AIC BIC

#> -62.2027 134.4054 145.4289

#>

#> Random effects covariance matrix:

#> StdDev

#> (Intercept) 0.1314838

#>

#> Fixed effects:

#> Estimate Std.Err z-value p-value

#> (Intercept) -0.3308 0.3136 -1.0547 0.291551

#> mhhx1 -0.9077 0.5221 -1.7387 0.082086

#> gender1 -0.3253 0.4960 -0.6559 0.511905

#> exposure1 -0.0532 0.4749 -0.1121 0.910769

#>

#> Integration:

#> method: adaptive Gauss-Hermite quadrature rule

#> quadrature points: 11

#>

#> Optimization:

#> method: hybrid EM and quasi-Newton

#> converged: TRUE

#Extract icc

icc(m)

#> # Intraclass Correlation Coefficient

#>

#> Adjusted ICC: 0.005

#> Conditional ICC: 0.005

#Simulate new data from model and fit new model to simulated data

nsim = 10

results = rep(NA, nsim)

for (isim in 1:nsim) {

newy=GLMMadaptive::simulate(m)

newfit = update(m, newy ~ .)

results[isim] = icc(newfit)[[1]]

}

#> Warning in mixed_model(fixed = distress ~ mhhx + gender + exposure, random = ~1

#> | : unknown names in control: formula

head(newy) # Visualise data

#> [,1]

#> [1,] 0

#> [2,] 1

#> [3,] 0

#> [4,] 0

#> [5,] 0

#> [6,] 0

str(newy)

#> num [1:100, 1] 0 1 0 0 0 0 0 0 1 0 ...

all(newy == m$data$distress) # Checking whether simulated data different to original data

#> [1] FALSE

summary(newfit) # Checking model (estimates are the same as m)

#>

#> Call:

#> mixed_model(fixed = distress ~ mhhx + gender + exposure, random = ~1 |

#> school, data = data, family = "binomial", formula = newy ~

#> mhhx + gender + exposure)

#>

#> Data Descriptives:

#> Number of Observations: 100

#> Number of Groups: 67

#>

#> Model:

#> family: binomial

#> link: logit

#>

#> Fit statistics:

#> log.Lik AIC BIC

#> -62.2027 134.4054 145.4289

#>

#> Random effects covariance matrix:

#> StdDev

#> (Intercept) 0.1314838

#>

#> Fixed effects:

#> Estimate Std.Err z-value p-value

#> (Intercept) -0.3308 0.3136 -1.0547 0.291551

#> mhhx1 -0.9077 0.5221 -1.7387 0.082086

#> gender1 -0.3253 0.4960 -0.6559 0.511905

#> exposure1 -0.0532 0.4749 -0.1121 0.910769

#>

#> Integration:

#> method: adaptive Gauss-Hermite quadrature rule

#> quadrature points: 11

#>

#> Optimization:

#> method: hybrid EM and quasi-Newton

#> converged: TRUE

#Extract paramaters from refitted model

icc(newfit)

#> # Intraclass Correlation Coefficient

#>

#> Adjusted ICC: 0.005

#> Conditional ICC: 0.005

quantile(results, c(0.025, 0.975))

#> 2.5% 97.5%

#> 0.005227451 0.005227451