*Bounty: 50*

*Bounty: 50*

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 2009^{1} (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:

gives`lmer`

*identical results for scaled and unscaled weights*, but different results with the expanded datasetgives`glmer`

*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

$$

logmathcal{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.}