Binomial GLM and different sample sizes

Question

I have a data set which consists of binomial proportions, let's say the success rate of converting a customer depending on the advertisement, the customer age, and various other factors.

For some common combinations of covariates, I have a lot of data, and therefore the binomial proportion of successes has low variance. For rare combinations of covariates, however, I have very little data, and therefore the variance of the proportion is high.

The magnitude of differences is very large, for example I might have 1 million trials for some combinations of covariates, and only 50 for others. However, I want to include ALL data in my model and weight it appropriately to get the best model fit.

I've tried to use R to do binomial (logistic) regression using a generalized linear model:

lrfit <- glm ( cbind(converted,not_converted) ~ advertisement + age, family = binomial)

This is a good start because it automatically weights the observations by the number of trials.

However, it's not good enough. Here's why: Let's say you have some observations with 100,000 trials and others with 1,000,000 trials. If you weight by number of trials the latter group is going to receive 10 times the weight. This seems nonsensical, however, because both observations are easily precise enough to receive equal treatment in the model. Clearly you want to penalize groups with only 10 or 100 trials, but as the number of trials gets larger, the weight should stop increasing.

Since in weighted least squares the reciprocal of the error variance is used as the weight, my idea would be to use calculate the posterior variance of the proportion (using Jeffrey's prior), then add some constant term to it (this will make sure the variance stops increasing at a certain number of trials) and then use the reciprocal of the sum as the weight.

Is this approach reasonable? Am I missing something? Can someone give me more information about this method?

Answer 1

Could you get the correct level of inferences from a Generalized Linear Mixed Model? Treating the advertisement as a practical random effect should introduce the correct level of penalization only on the smaller sampled levels. Those levels with any reasonable magnitude will get the appropriate fixed effect estimates.

Assuming the use of lme4 in R, I'm not sure if you'd need to replicate the rows or not. (I'm leaning towards yes). Then you would have something like:

glmm.fit <- glmer(success~ns(age,df=6)+(1|advertisement),family=binomial(),data=rep.df)

(Notice I splined the age using ns() from splines package since I can't imagine it would actually be linear.)

Answer 2

You may want to consider the arcsine transformation of your response.

If

$X \sim$ Bin($n,p$)

then

$E(\sin^{-1} (\sqrt{\frac{X+c}{n+2c}})) = \sin^{-1} \sqrt{p} + p^{-1/2}O(c-1/2) + O(p^{-3/2})$

and

$n Var(sin^{-1}(\sqrt{\frac{X+c}{n+2c}}) = 1/4 + p^{-1}O(c-3/8) + O(p^{-2})$

Set $c=3/8$ and you have the variance stabilizing transformation. Then you could use a weighted least squares model where $n$ is your weight under the assumption that your variance covariance matrix is heteroskedastic, but diagonal.

Answer 3

With so many observations I would be worried about extra-binomial variance, that is, overdispersion. Did you try with a quasibinomial family?

That is, change your R code to

lrfit.quasi <- glm ( cbind(converted,not_converted) ~ advertisement + age, family = quasibinomial)

Can you tell us what happens when you do that? What was the estimated overdispersion factor? I would guess it could be much larger than one, the binomial case. Also, with 1 million observations, I guess there could be some extra structure to the sampling, like different sampling sites, subgroups, differing structure with time, autocorrelation ... That would introduce dependencies which would have to be modeled. You need to tell us much more about how you obtained this data.

The binomial model also assumes that the count is exact, with 1 million counts that could easily be false ... That is another source of extrabinomial variance! An alternative to using a quasibinomial model is to use some mixed model (or a negative binomial model).

Say, with $n=1000000$ and binomial $p$ around 0.5, then the (binomial) standard error of $\hat{p}$ is around $0.0005$. An error in the count of about $500$ then changes the estimated $p$ with about the same amount, one standard error. That is a count error of about $1 ^0\!\!/\!_{00}$, one in thousand. Even a quite accurate count could give such and larger errors! Are you really sure your counts are without error?

If you would investigate what would happen with the fit if you try to model errors in the counts, the last chapter (on optimization) in Venables & Ripley MASS (fourth edition) have some theory and methods for that case.