3
$\begingroup$

was tasked with developing a regression model looking at student enrollment in different programs. This is a very nice, clean data set where the enrollment counts follow a Poisson distribution well. I fit a model in R (using both GLM and Zero Inflated Poisson.) The resulting residuals seemed reasonable.

However, I was then instructed to change the count of students to a "rate" which was calculated as students / school_population (Each school has its own population.)) This is now no longer a count variable, but a proportion between 0 and 1. This is considered the "proportion of enrollment" in a program.

This "rate" (students/population) is no longer Poisson, but is certainly not normal either. So, I'm a bit lost as to the appropriate distribution, and subsequent model to represent it.

A log normal distribution seems to fit this rate parameter well, however I have many 0 values, so it won't actually fit.

Any suggestions on the best form of distribution for this new parameter, and how to model it in R?

$\endgroup$
2

3 Answers 3

4
$\begingroup$

Most software that supports Poisson regression will support an offset and the resulting estimates will become log(rate) or more acccurately in this case log(proportions) if the offset is constructed properly:

 # The R form for estimating proportions
 propfit <- glm( DV ~ IVs + offset(log(class_size), data=dat, family="poisson")

exp(coef(propfit) : should be the baseline proportion exp(Intercept) followed by the factors to multiply that baseline proportion estimate to arrive at the class-specific estimates.

$\endgroup$
3
  • $\begingroup$ that is the solution I've seen elsewhere. For some reason, that regression doesn't seem to fit as well when I include the offset. In fact, the regression fails for many of the IVs. (I'm iterating over different IVs to find the best ones in a form of stepwise regression.) $\endgroup$
    – Noah
    Commented Apr 17, 2013 at 21:36
  • $\begingroup$ One of the issues is how to use/interpret the predict.glm() funciton. Without an offset, I can just use predict.glm(fit, type="response") and get the predicted values. With an offset, how does the prediction differ? The help page for predict.glm() isn't clear on this. $\endgroup$
    – Noah
    Commented Apr 17, 2013 at 21:37
  • $\begingroup$ It will be giving you the predicted proportions and you should check this by multiplying the N by the predicted values. $\endgroup$
    – DWin
    Commented Apr 17, 2013 at 21:42
2
$\begingroup$

If the rate is a count in a given category divided by so many in total, you could do a binomial GLM; there you're in effect explicitly modelling a proportion...

glm( formula, family=binomial)

The default link is logit rather than log as for the Poisson, but for small proportions the difference is negligible.

If you're looking at many categories, there are a number of extensions that might be appropriate, such as multinomial logit choice models.

Also see http://en.wikipedia.org/wiki/Generalized_linear_model#Multinomial_regression

$\endgroup$
3
  • $\begingroup$ A rate is never a count. $\endgroup$
    – DWin
    Commented Apr 17, 2013 at 21:42
  • 1
    $\begingroup$ @DWin I have replaced the phrase "out of" by "divided by" in order to make the original intent clearer. A 'count divided by' something can be a rate. The original 'out of' was intended to convey the same meaning while emphasizing the connection between rates, and proportions-in-the-binomial (but was, I admit, unclear to the point of being potentially misleading). $\endgroup$
    – Glen_b
    Commented Apr 17, 2013 at 23:09
  • $\begingroup$ A rate is a count divided by an exposure which is in turn a sum of time*counts; a proportion is a count divided by a (larger) count of individuals at risk. A logit is a logged( count divided by a "non-count"). $\endgroup$
    – DWin
    Commented Oct 20, 2014 at 22:57
0
$\begingroup$

I don't get why you think you couldn't transform that DV to log-normality only because of the zeros. You can just add a constant to this variable and then log-tranform it. Usual numbers to do so are 1, 10 or the minimum value of your variable. It would therefore have the form of log(1+DV), log(10+DV), log(min(DV)+DV) or in a more general term, log(a+DV). In R you can use the log1p() function instead of log() to automatically apply a log(1+DV).

If it doesn't change the distribution of he residuals in your models, you can also try other transformations, such as DV^a (in which a represents any number, usually 2, 3 or 0.5) or even exp(DV).

Also, another interesting option is to consider using robust regression, which can handle bad behaved distribugions.

For additional help, bringing us a densityplot() image of you residuals would be great.

$\endgroup$
1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.