I have been looking into how to use the weight argument of glmer
/lmer
to represent "frequency" weights. I was comparing model estimations with an expanded dataset, with the frequency weights inputted directly to the weight argument, and with the weights scaled as in Method A in Carle 20091 (i.e. scaled such that their sum matches the number of observations). I have found a puzzling difference in how they affect the results.
Specifically:
lmer
gives identical results for scaled and unscaled weights, but different results with the expanded datasetglmer
gives identical results with the unscaled weights as with the expanded dataset, but different results with the scaled weights
Could anybody clarify what is the reason for the difference?
Backgroud
I want to compare performance of two groups on a cognitive psychology task. The groups have been matched by a collaborator according to relevant covariates (e.g. socioeconomic status) and as a result few participants from the control group have a frequency weight of 2, indicating that they should be counted twice in the analyses in order for the groups to be matched on all the covariates. I would like to use multilevel models, and I am trying to see if I can use the frequency weight with the weight argument of lme4
(Note: I am not sure of whether frequency weight is the correct term here; perhaps sampling weight is more appropriate). The dependent variable is binomial, so we'd need to use a GLMM.
Toy example:
Generate some data
# covariance matrix
sigma = matrix(c(1, 0.7, 0.7, 1), ncol = 2)
# simulate data
set.seed(4321)
df = mvrnorm(n = 100, mu = c(0,0), Sigma = sigma)
# Simulate clustering and calculate weights
df = data.frame(df)
names(df) = c("Y", "X")
df$cluster = rep(1:10, each = 10)
df = within(df, {
Y = Y + cluster/10
X = X + cluster/10
weights = rep(1:2, each = 50)
scaled_weights = weights*2/3
})
#data set with double weighted observations entered twice
df2 = rbind(df, df[51:100,])
Effects of weights in lmer
. Here the models with weights gives the same results, regardless of the scaling, and the model with the expanded dataset gives different estimates.
# Models
m.weights = lmer(Y ~ X + (1 | cluster), df, weights = weights) # weights
m.double = lmer(Y ~ X + (1 | cluster), df2) # expanded dataset
m.scaled = lmer(Y ~ X + (1 | cluster), df, weights = scaled_weights) #scaled weights
round(summary(m.weights)$coef,digits=4)
# Estimate Std. Error t value
# (Intercept) 0.2498 0.0904 2.7645
# X 0.6908 0.0658 10.4928
round(summary(m.scaled)$coef,digits=4) # SAME AS m.weights
# Estimate Std. Error t value
# (Intercept) 0.2498 0.0904 2.7645
# X 0.6908 0.0658 10.4928
round(summary(m.double)$coef,digits=4) # DIFFERENT
# Estimate Std. Error t value
# (Intercept) 0.2445 0.0842 2.9029
# X 0.6847 0.0538 12.7332
Effects of weights in glmer
# discretize dependent variable in 2 levels
df$dY <- ifelse(df$Y>0.8,1,0)
df2$dY <- ifelse(df2$Y>0.8,1,0)
m2.weights = glmer(dY ~ X + (1 | cluster), df, weights = weights, family =binomial) # weights (unscaled)
m2.double = glmer(dY ~ X + (1 | cluster), df2, family =binomial) # expanded dataset
m2.scaled = glmer(dY ~ X + (1 | cluster), df, weights = scaled_weights, family =binomial) # scaled weights
# Warning message:
# In eval(family$initialize, rho) : non-integer #successes in a binomial glm!
Note that the model with scaled weights give a warning, because in glm(er) the weights should correspond to the number of trials when the dependent variable is the proportion of successes. However as I understand them the weights are taken into account in the log-likelihood as
$$ \log\mathcal{L}(\boldsymbol{\theta}) = \sum_{i = 1}^{n} w_i \log \left[ p(y_i | \boldsymbol{x}_i,\boldsymbol{\theta})\right] $$
where $w_i$ is the weight. If that's correct, each observation $i$ contributes $w_i$ times its log-likelihood to the total log-likelihood, which should be the same as saying that each $i$ represents $w_i$ observations. (Since this is what I would like to express with the weights, I thought I could use this argument for my purposes). Furthermore, based on this I would expect to find that the model with the unscaled weight and that with the expanded dataset to give similar or identical results, and indeed that is what I find - in this case is the model with scaled weights that differs from the other two.
round(summary(m2.weights)$coef,digits=4)
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -1.3028 0.3902 -3.3390 8e-04
# X 1.5741 0.2955 5.3263 0e+00
round(summary(m2.scaled)$coef,digits=4) # DIFFERENT
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -1.2620 0.4183 -3.0167 0.0026
# X 1.5579 0.3531 4.4114 0.0000
round(summary(m2.double)$coef,digits=4) # SAME as m.weights
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -1.3028 0.3902 -3.3390 8e-04
# X 1.5741 0.2955 5.3263 0e+00
1 Carle, A. C. (2009). Fitting multilevel models in complex survey data with design weights: Recommendations. BMC Medical Research Methodology, 9 (1), 1–13.