Random-effects and proportions

Question

I'm having trouble developing an approach to modeling the longitudinal effect of an intevention on changes in vaccination rates.

The data come from an annual collection in 600 villages of information regarding the number of individuals aged 12 years in each village and each individual's dichotomous (vaccination) status. So, for each village, there is an annual count of positive status and negative status 12 year olds. The counts are available for 7 years. Individuals were not recounted. They are effectively replaced each year in each village by a new cohort of children aged 12. After the third year and before the year 4 count, an intervention was introduced in a random selection of 300 villages.

The number of 12 year old children in each village varies from year; there is also considerable variation in the size of the villages.

The main question of interest is whether the proportion of 12 year old children with positive vaccination status is increasing/descreasing/stable as a consequence of the intervention.

I have been singularly unsuccessful in finding any similar analysis, or a tutorial on such analysis, nor even a similarly constructed set of data. So I'd appreciate any guidance related to that, or to a suitable text to read.

Answer 1

This is easily done using a generalized linear mixed model. I'll be using R code so I hope you're familiar with that. If not, I'm sure all the other major software packages can handle this.

You have a repeated measures design (panel data) so within a regression model you need a random effect to take this into account. Village will be the random effect, so this takes into account the specific base probability of vaccination within each village (a village that has high proportions of vaccinated 12-year olds will probably have high proportions the next year too):

(1|village)

So this means we have to use a mixed model. The fixed effects might be intervention and time (perhaps there is a trend over time?) and perhaps also their interaction term:

Intervention * Time

The dependent variable is a proportion, and if you have the number of 12 year-olds in each village for each measurement occasion (and we can assume no overlap), we can use a binomial distribution function by using the weights command:

library(lmer)
glmer(proportion ~ intervention * time + (1|village), weights=size, family=binomial)

And this is it. To my understanding the regression model will treat each village as a number of rows equal to the size variable (the number of 12 year-olds in the village at that time), were some score 0 (no vaccination) and some score (vaccination), and the proportion is equal to the proportion variable in your data.

Say that the first village at the first measurement has 300 12 year-olds, and 20% are vaccinated (proportion = 0.2). Using the above syntax is equivalent to having a data set with 300 rows for that village and measurement point, where 60 of the 300 rows have a '1' for vaccination and the other 240 have a '0'.

I hope this helps.

EDIT:

Another option is to use a general linear mixed model with Poisson distribution. I haven't seen this used in a mixed model, but I can't see any reason why it shouldn't work, and I did some simulation experiments that indicate that the results were as expected.

The idea is that poisson regression can be used to model count data, and if we want to model a proportion, we can express this as a count out of a population. Say you have a village with 100 12-year olds, and 23% have gotten the vaccination. In the example above, you the proportion variable would have the value 0.23 for this row, and the size variable would be 100. This is equivalent to saying that there was a count of 23 vaccinated individuals out of 100. So if you multiply proportion with size for each row, you get the counts.

There is an excellent post here, by the user ocram, that explains what happens next. I'll just copy and paste his text here since it explains the process so succinctly and I can't explain it as well as this:

Poisson regression is typically used to model count data. But, sometimes, it is more relevant to model rates instead of counts. This is relevant when, e.g., individuals are not followed the same amount of time. For example, six cases over 1 year should not amount to the same as six cases over 10 years. So, instead of having

$\log \mu_x = \beta_0 + \beta_1 x$

(where $\mu_x$ is the expected count for those with covariate $x$), you have

$\log \tfrac{\mu_x}{t_x} = \beta'_0 + \beta'_1 x$

(where $t_x$ is the exposure time for those with covariate $x$). Now, the last equation could be rewritten

$\log \mu_x = \log t_x + \beta'_0 + \beta'_1 x$ and $\log t_x$ plays the role of an offset.

In this context, it's applied to individuals followed over different lengths of time, but it applies to your situation as well. You thus want to model the logarithm of events divided by size, and you accomplish this by adding an offset variable of log(size) on the right side of the equation. In R notation:

vaccinated.count <- proportion * size
glmer(vaccinated.count ~ offset(log(size)) + intervention * time + (1|village), family=poisson)

The estimate for intervention will indicate the (log) change in proportion of vaccinated 12-year olds in villages which have gotten the intervention rather than the odds ratio in the model based on the binomial distribution.