# Binomial GLM and different sample sizes

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?

• I don't understand this line "both observations are easily precise enough to receive equal treatment in the model": one has 10 times the observations than the other, and so has a 1/10th of the variance. Granted these variances will be small, but one is still smaller than the other. – Simon Byrne Mar 15 '13 at 16:34
• Yes, that wasn't very clear. What I meant was that although the two observations have different variances, it intuitively doesn't seem right to give one of them 10 times the weight of the other, if you consider the fact that the model is not going to be accurate enough for such small differences to matter anyway. Adding a billion more trials when you're already very confident about a proportion simply shouldn't result in major parameter changes. – M. Cypher Mar 15 '13 at 16:40
• Does down-weighting actually result in major parameter changes? By the way if you want to re-weight, you can just use the weights argument in glm, which will automatically deal with the fact that different probabilities give different variances. – Simon Byrne Mar 16 '13 at 13:47
• You can use any weights you like, but you are wrong that the standard approach is "nonsensical." Even a huge ratio of weights will not entirely suppress the influence of the low-weighted cases, especially if they have unique or unusual combinations of independent variable values. If introducing more data creates "major parameter changes," then you have a more fundamental problem: either the model is a bad fit and/or it changes over time. Fiddling with weights won't address that problem; it will only give you false confidence in the results. – whuber Jun 13 '13 at 18:17
• You seem to be worried about the weights. Are you able to produce some data which confirms your theory? For example, choose a "true beta", choose some covariates, simulate binomial data using these probabilities (using n's ranging from $5$ up to $1000000$). Then fit the glm and show it is producing "non sensical" results. – probabilityislogic Jan 12 '14 at 11:44

## 3 Answers

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.)

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.

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.

• I think using the arcsine transformation would be going backwards since the GLM with logit link is more accurate and there's no need to stabilize variances. I can set the weights in my model to anything I like, what I'm looking for is an explanation of what weighting would be the optimal one, and why. – M. Cypher Mar 16 '13 at 16:06
• Using the arcsine trenaformation (followed by linera least squares fitting) can be useful, especially in simple models. It will also work correctly with overdispersed data! – kjetil b halvorsen Apr 24 '16 at 20:29