Inspect the estimates and calculate \(R^2_m\)

The summ function from the package jtools (Long, 2020) gives visually pleasing summaries of fitted lme4 models and reports the marginal effect size for the fixed effects (\(R^2_m\)). Additionally, this function reports p-values based on the Satterthwaite approximation for the degrees of freedom. The lme4 package deliberately omits reporting p-values due to issues regarding the calculation of the degrees of freedom (for details on evaluating the fixed effects by using Null Hypothesis Significance Testing approaches, see, Luke, 2017).

summ(model.1) 

The fixed effect for standlrt is estimated to be larger than zero, which is evident both by its SE and by the highly significant p-value (rejecting the null hypothesis which states that the fixed effect is equal to zero). The value for the ICC in this data set is .14, which indicates the amount of within-school clustering with respect to normexam, before accounting for (part of) the explained variance by introducing (random) predictors.

fixef <- fixef(model.1)  # extract the fixed effect 
y_hat <- fixef[1] + fixef[2]*tutorial$standlrt # predict the outcome 

sigma_f  <- var(y_hat) # obtain the variance in the predicted outcome based on the fixed effect

# extract the estimated random effects from the fitted model (you can use the summary function)
# i.e., summary(model.1)
sigma_u0 <- 0.09044
sigma_u1 <- 0.01454
sigma_e  <- 0.55366

# calculate the R^2_m
marginal_Rsq_fixed <- (sigma_f)/(sigma_f + sigma_u0 + sigma_u1 + sigma_e)
marginal_Rsq_fixed

## [1] 0.3170484

Compare the models

Since we are using uninformative priors, we will take a small detour and show that the Bayesian posterior parameter estimates correspond to those obtained with lmer when using FML.

# Check and inspect the model
model.def <- jags.model(file = textConnection(jags.model),
                        inits = list(.RNG.name="base::Wichmann-Hill",
                                     .RNG.seed=100),
                        data = tutorial, n.chains = 2)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 4059
##    Unobserved stochastic nodes: 72
##    Total graph size: 16539
## 
## Initializing model

update(object = model.def, n.iter = 1000)

# ask only for these parameters
parameters <- c("mu.alpha", "mu.beta", "tau", "invSigma")

results <- coda.samples(model = model.def, variable.names = parameters, n.iter =1000)
summary(results)

## 
## Iterations = 2001:3000
## Thinning interval = 1 
## Number of chains = 2 
## Sample size per chain = 1000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                    Mean       SD  Naive SE Time-series SE
## invSigma[1,1]  13.52115  3.58318 0.0801224      0.1558283
## invSigma[2,1] -16.27213  8.90956 0.1992239      0.5807041
## invSigma[1,2] -16.27213  8.90956 0.1992239      0.5807041
## invSigma[2,2]  79.44014 31.90626 0.7134457      2.1708230
## mu.alpha       -0.01271  0.04361 0.0009752      0.0013995
## mu.beta         0.55412  0.02121 0.0004743      0.0008501
## tau             1.80630  0.04071 0.0009103      0.0009105
## 
## 2. Quantiles for each variable:
## 
##                    2.5%       25%       50%       75%     97.5%
## invSigma[1,1]   8.05983  11.05483  13.07467  15.39220  22.17893
## invSigma[2,1] -39.38018 -20.46357 -14.63042 -10.09460  -4.02706
## invSigma[1,2] -39.38018 -20.46357 -14.63042 -10.09460  -4.02706
## invSigma[2,2]  37.66881  57.83871  72.84995  93.21918 170.70371
## mu.alpha       -0.09981  -0.04166  -0.01324   0.01561   0.06769
## mu.beta         0.51329   0.53977   0.55368   0.56767   0.59710
## tau             1.72764   1.77944   1.80671   1.83268   1.89180

#transform the precisions of the random eff. into standard deviations 
sqrt(1/13.52115) # intercept SD

## [1] 0.2719526

sqrt(1/79.44014) # slope SD

## [1] 0.1121967

sqrt(1/1.8063)   # residual SD

## [1] 0.744055

Note how the posterior mean estimates along with their standard deviations roughly correspond to their (Full) Maximum Likelihood counterparts. Note, the (co)variance components are expressed in terms of precisions, the reason for this are beyond the scope of the tutorial. Thus, in order to be able to compare the results with the ones from lmer, we need to take the inverse for the values of invSigma and tau, and then take the square root (since the random effects summarized by the summ function are expressed as standard deviations). The last three rows from the output contain the estimated random effects expressed in terms of standard deviations.

Obtain multiple imputed data sets

Now, we estimate the model again and save all the sampled random effects, which we will treat as multiple imputed values. We run the sampler by using two chains with a 1000 iterations each (excluding the burn-in period of a 1000 extra iterations).

# fit the model again but ask to monitor all the random effects
model.def <- jags.model(file = textConnection(jags.model),
                        inits = list(.RNG.name="base::Wichmann-Hill",
                                     .RNG.seed=100),
                        data = tutorial, n.chains = 2)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 4059
##    Unobserved stochastic nodes: 72
##    Total graph size: 16539
## 
## Initializing model

update(object = model.def, n.iter = 1000)

# the random effects (the ones we need for the aim of this approach)
parameters <- c("alpha", "beta", "mu.alpha", "mu.beta")

results <- coda.samples(model = model.def, variable.names = parameters, n.iter =1000)

Extract the samples from the posterior (from both chains):

chain1 <- as.data.frame(results[[1]])
chain2 <- as.data.frame(results[[2]])
samples <- rbind(chain1, chain2) # combine both chains

Before we start manipulating the sampled random (and fixed) effects, such that we can add each one to a copy of the original data set, it would be nice to have a glimpse of the data set’s layout.

datatable(tutorial.1, options = list(pageLength = 5))

Extract the fixed effect from each sampled vector (i.e., \(\alpha\) and \(\beta)\):

fixed_alphas <- samples[, 131]
fixed_betas <- samples[, 132]
fixed_alphas <- as.data.frame(fixed_alphas)
fixed_betas <- as.data.frame(fixed_betas)

Repeat each fixed effect N times (in this case 4059 - the number of \(N_{level-1}\) observations) and store it in a separate data frame that will be merged with the imputed data sets later on:

# the fixed intercepts
samp_fixed_alphas <- list()
for (i in 1:nrow(fixed_alphas)){
   samp_fixed_alphas[[i]] <-rep.int(fixed_alphas[i, 1], nrow(tutorial.1))  
}

# the fixed slopes
samp_fixed_betas <- list()
for (i in 1:nrow(fixed_betas)){
   samp_fixed_betas[[i]] <-rep.int(fixed_betas[i, 1], nrow(tutorial.1))  
}

Get the sampled random intercepts (\(\alpha_j's\))

First, split the data per group:

split <- split(tutorial.1, tutorial.1$school)

Afterwards, extract the size of each group:

n_gr <- sapply(split, nrow) 

Exclude the columns containing the fixed effects from the samples data set:

samples <- samples[, -c(131,132)]

Extract the random \(\alpha_{j}s\) and \(\beta_{j}s\) into separate data frames:

alphas <- samples [, 1:65]  # extract the random intercepts
betas <- samples [, 66:130] # extract the random slopes

Use the size of each group to replicate, from each iteration, the respective \(\alpha_j\), n_gr number of times:

samp_alphas <- list()

for (i in 1:nrow(alphas)){
   samp_alphas[[i]] <-rep.int(alphas[i, ], n_gr)  
}

Extract these in a big matrix with the samples from every iteration stored in a separate column:

samp_alphas_mat <- matrix(nrow = nrow(tutorial.1), ncol = nrow(alphas))

for(i in 1:nrow(alphas)){
  samp_alphas_mat[, i] <- unlist(samp_alphas[[i]])
}

Get the sampled random slopes (\(\beta_j's\))

Use the size of each group to replicate, from each iteration, the respective \(\beta_j\), n_gr number of times:

samp_betas <- list()

for (i in 1:nrow(betas)){
   samp_betas[[i]] <-rep.int(betas[i, ], n_gr)  
}

Extract them in a big matrix with the samples from every iteration stored in a separate column:

samp_betas_mat <- matrix(nrow = nrow(tutorial.1), ncol = nrow(betas))

for(i in 1:nrow(betas)){
  samp_betas_mat[, i] <- unlist(samp_betas[[i]])
}

Combine the sampled random effects

First, combine the \(\alpha_j's\):

imputed <- list()

for (i in 1:ncol(samp_alphas_mat)){
  imputed[[i]] <- cbind(tutorial.1, samp_alphas_mat[, i])
}

Add the \(\beta_j's\):

for (i in 1:ncol(samp_betas_mat)){
  imputed[[i]] <- cbind(imputed[[i]], samp_betas_mat[, i])
}

Add the sampled fixed effects:

The \(\alpha's\):

for (i in 1:length(imputed)){
  imputed[[i]] <- cbind(imputed[[i]], samp_fixed_alphas[[i]])
}

The \(\beta's\):

for (i in 1:length(imputed)){
  imputed[[i]] <- cbind(imputed[[i]], samp_fixed_betas[[i]])
}

Voilà, we obtained 2000 imputed data sets.

Transform the outcome variable

In order to be able to fit multiple linear regression models we need to transform the outcome variable as follows,

\[Z_{ij} = \text{normexam}_{ij} - \alpha_j - \beta_j\;\text{standlrt}_{ij} + \alpha + \beta\;\text{standlrt}_{ij},\]

which leads to \[ Z_i = \eta_0 + \eta_1\;\text{standlrt}_i + e_i.\]

Construct a function that calculates \(Z\):

transform_z <- function(df){
  df$z <- df[, 3] - df[, 5] - df[, 6] * df[, 4] + df[, 7] + df[, 8] * df[, 4]
  df
}

Apply the function to the “imputed” data sets and obtain a column for \(Z\):

imputed <- lapply(imputed, transform_z)

Now let’s quickly have a look at one of the imputed data sets:

datatable(imputed[[1]][,-c(1,2)], options = list(pageLength = 5)) # exclude the school and student columns 

Estimates:

Fit linear regression models to each imputed data set and extract the estimated coefficients and their respective variance-covariance matrices:

estimates <- list()
vCOV <- list()
for(i in seq_along(imputed)){
  estimates[[i]] <- coef(lm(imputed[[i]][,9] ~ imputed[[i]][,4]))
  vCOV [[i]] <- vcov(lm(imputed[[i]][,9] ~ imputed[[i]][,4]))
}

Extract these estimates in a data frame:

estimates <- unlist(estimates)
intercepts <- estimates[seq(1,length(estimates),2)] # select every other element
slopes <- estimates[seq(2,length(estimates),2)]     # opposite
estimates <- data.frame(intercepts, slopes)         # combine

Apply the multiple imputation equations

The following equations are taken from Van Buuren (2018, Ch. 2.3).

m <- nrow(estimates) # number of "imputations" (i.e., iterations in this case)

Combined estimate:

\[ \bar{\boldsymbol{\eta}} = \frac{1}{m} \sum_{l = 1}^{m} \hat{\boldsymbol{\eta}_l} \]

eta_bar <- apply(estimates, 2, mean) 
eta_bar <- t(eta_bar)
eta_bar <- as.matrix(eta_bar)

Average over the variances:

\[ \bar{U} = \frac{1}{m} \sum_{l = 1}^{m} \bar{U}_l \]

Since the list vCOV contains variance-covariance matrices and we need to take the average across all of them (and not within them), it is best to extract every element of each covariance matrix and then take the average:

V <- matrix(nrow = nrow(estimates), ncol = 4)

for(i in seq_along(vCOV)){
  V[i, 1] <- vCOV[[i]][1,1]
  V[i, 2] <- vCOV[[i]][1,2]
  V[i, 3] <- vCOV[[i]][2,1]
  V[i, 4] <- vCOV[[i]][2,2]
}

U_bar <- apply(V, 2, mean)
U_bar <- matrix(U_bar, nrow = 2, ncol = 2)

Unbiased estimate of the variance between the m complete data estimates:

\[B = \frac{1}{m-1} \sum_{l = 1}^{m} (\hat{\boldsymbol{\eta}}_l - \bar{\boldsymbol{\eta}})(\hat{\boldsymbol{\eta}}_l - \bar{\boldsymbol{\eta}})'.\]

B <- cov(estimates)

Total variance:

\[T = \bar{U} + (1 + \frac{1}{m})*B\]

Total_var <- U_bar + (1 + 1/m)*B

Proportion of variation attributable to the missing data (a compromise over all estimates)

\[\bar{\lambda} =(1 + \frac{1}{m}) tr(BT^{-1})/k\]

k <- ncol(estimates) # number of parameters in eta_hat
lambda_hat <- (1+ 1/m) * sum(diag(B %*% ginv(Total_var)))/k

Degrees of freedom:

\[\nu_{old} = \frac{m-1}{\bar{\lambda}^2}; \; \nu_{com} = N_{level-1} - k; \; \nu_{obs} = \frac{\nu_{com} + 1}{\nu_{com} + 3}\nu_{com}(1-\bar{\lambda}); \; \nu = \frac{\nu_{old} \nu_{obs}}{\nu_{old} + \nu_{obs}}.\]

N <- nrow(tutorial.1) # N_level-1

nu_old <- (m - 1)/lambda_hat^2

nu_com <- N - k

nu_obs <- ((nu_com +1)/(nu_com + 3))*nu_com*(1 - lambda_hat)

nu <- (nu_old*nu_obs)/(nu_old + nu_obs)

Fraction of missing information:

\[\gamma = \frac{\nu + 1}{\nu + 3}\bar{\lambda}+\frac{2}{\nu + 3}\]

gamma <- ((nu + 1)/(nu + 3)) * lambda_hat + (2/(nu + 3))

Calculate the MI-based \(N_{eff}\)

\[\text{MI-based}\; N_{eff} = N_{level-1} - \gamma*N_{level-1}\]

MI_N_eff <- N - gamma*N  
MI_N_eff

## [1] 874.2342

The MI-based effective sample size is 874.

Compare with the ICC-based \(N_{eff}\)

Calculate the ICC-based \(N_{eff}\)

We already know that the ICC = .17 since it was given by the summ function, however, here we calculate it by hand, to illustrate that the ICC can only be computed using the estimated variances from the intercept-only model and the average school (group) size:

\[ICC = \frac{\sigma_{u0}}{\sigma_{u0} + \sigma_{\epsilon}}, \] where \(\sigma^2_{u0}\) denotes the variance for the random slope and \(\sigma^2_{\epsilon}\) denotes the residual variance estimated from a random intercept-only model.

# fit a random intercept-only model
model.0 <- lmer(normexam ~1 + (1|school), REML = FALSE, data = tutorial.1) 
summary(model.0)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: normexam ~ 1 + (1 | school)
##    Data: tutorial.1
## 
##      AIC      BIC   logLik deviance df.resid 
##  11016.6  11035.6  -5505.3  11010.6     4056 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9471 -0.6486  0.0117  0.6992  3.6578 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  school   (Intercept) 0.1686   0.4107  
##  Residual             0.8478   0.9207  
## Number of obs: 4059, groups:  school, 65
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) -0.01317    0.05363  -0.246

ICC <- 0.1686 / (0.1686 + 0.8478)
ICC

## [1] 0.1658796

The ICC-based \(N_{eff}\) is defined as

\[\text{ICC-based} \; N_{eff} = \frac{N_{level-1}}{1+(n_c - 1)ICC},\]

where \(n_c\) denotes the within-group sample size. As already mentioned, a drawback of the ICC-based \(N_{eff}\) is that when the data has varying group sizes, a compromise, such as taking the average group size, has to be made.

n_clus <- mean(sapply(split, nrow)) # the average group (i.e., school) size (62.4 in this case)
ICC_N_eff <- N / (1 + (n_clus - 1) * ICC)
ICC_N_eff 

## [1] 362.6483

The value of the ICC-based effective sample size is 363. This demonstrates that by accounting for the within-group variation through a model that includes a random intercept and random slope for the predictor standlrt, the number of effective level-1 observations increases from 362 to 874. Essentially, the model with the predictor “explains away” some of the clustering within groups indicated by the ICC. As an exercise, you may want to try recalculating the MI-based effective sample size by using a random intercept-only model with JAGS, and compare the value obtained to that of the ICC-based approach.

Inspect the parameter estimates and calculate \(R^2_m\)

summ(model.2) 

Note that the value for the fixed effect of standlrt is estimated to be zero (this can also be seen from the highly non-significant p-value). Moreover, now \(R^2_m = 0\), and we expect that our new default BF will tend towards 19.

Additionally, the ICC is .15 which yields an ICC based \(N_{eff}\) of 375 and MI-based \(N_{eff}\) of 1070.