10
$\begingroup$

In R, there are three methods to format the input data for a logistic regression using the glm function:

  1. Data can be in a "binary" format for each observation (e.g., y = 0 or 1 for each observation);
  2. Data can be in the "Wilkinson-Rogers" format (e.g., y = cbind(success, failure)) with each row representing one treatment; or
  3. Data can be in a weighted format for each observation (e.g., y = 0.3, weights = 10).

All three approach produces the same coefficient estimates, but differ in the degrees of freedom and resulting deviance values and AIC scores. The last two methods have fewer observations (and therefore degrees of freedom) because they use each treatment for the number of observations whereas the first uses each observation for the number of observations.

My question: Are there numerical or statistical advantages to using one input format over another? The only advantage I see is not having to reformat one's data in R to use with the model.

I have looked at the glm documentation, searched on the web, and this site and found one tangentially related post, but no guidance on this topic.

Here is a simulated example that demonstrates this behavior:

# Write function to help simulate data
drc4 <- function(x, b =1.0, c = 0, d = 1, e = 0){
    (d - c)/ (1 + exp(-b * (log(x)  - log(e))))
}
# simulate long form of dataset
nReps = 20
dfLong <- data.frame(dose = rep(seq(0, 10, by = 2), each = nReps))
dfLong$mortality <-rbinom(n = dim(dfLong)[1], size = 1,
                              prob = drc4(dfLong$dose, b = 2, e = 5))

# aggregate to create short form of dataset
dfShort <- aggregate(dfLong$mortality, by = list(dfLong$dose), 
                     FUN = sum)
colnames(dfShort) <- c("dose", "mortality")
dfShort$survival <- nReps - dfShort$mortality 
dfShort$nReps <- nReps
dfShort$mortalityP <- dfShort$mortality / dfShort$nReps

fitShort <- glm( cbind(mortality, survival) ~ dose, 
                 data = dfShort, 
                 family = "binomial")
summary(fitShort)

fitShortP <- glm( mortalityP ~ dose, data = dfShort, 
                  weights = nReps,     
                  family = "binomial")
summary(fitShortP)

fitLong <- glm( mortality ~ dose, data = dfLong, 
                family = "binomial")
summary(fitLong)
$\endgroup$
  • 1
    $\begingroup$ In your example, the difference between the null and residual deviance is equal for all three models. If you add or remove a parameter, the change in AIC is also the same for all three. $\endgroup$ – Jonny Lomond Jan 12 '18 at 12:53
  • 1
    $\begingroup$ You answered yourself: if using the same data, to same results, then how could they be different? Moreover, if the method would give different results only because of providing the data in different format, something would be deeply wrong with it, or with its implementation. $\endgroup$ – Tim Jan 12 '18 at 14:55
  • $\begingroup$ The WR format is ultimately a weighted likelihood. The problem with weights is R can't tell whether they're frequency weights, probability weights, or other. With survey weighting, e.g., you may only have n observations but they represent segments of the population/sampling frame. So the degrees of freedom is indeed 100. svyglm from the survey package gives you better methods of handling the weight argument. $\endgroup$ – AdamO Jan 12 '18 at 19:04
7
+50
$\begingroup$

There's no statistical reason to prefer one to the other, besides conceptual clarity. Although the reported deviance values are different, these differences are completely due to the saturated model. So any comparison using relative deviance between models is unaffected, since the saturated model log-likelihood cancels.

I think it's useful to go through the explicit deviance calculation.

The deviance of a model is 2*(LL(Saturated Model) - LL(Model)). Suppose you have $i$ different cells, where $n_i$ is the number of observations in cell $i$, $p_i$ is the model prediction for all observations in cell $i$, and $y_{ij}$ is the observed value (0 or 1) for the $j$-th observation in cell $i$.

Long Form

The log likelihood of the (proposed or null) model is $$\sum_i\sum_j\left(\log(p_i)y_{ij} + \log(1 - p_i)(1 - y_{ij})\right)$$

and the log likelihood of the saturated model is $$\sum_i\sum_j \left(\log(y_{ij})y_{ij} + \log(1 - y_{ij})(1-y_{ij})\right).$$ This is equal to 0, because $y_{ij}$ is either 0 or 1. Note $\log(0)$ is undefined, but for convenience please read $0\log(0)$ as shorthand for $\lim_{x \to 0^+}x\log(x)$, which is 0.

Short form (weighted)

Note that a binomial distribution can't actually take non-integer values, but we can nonetheless calculate a "log likelihood" by using the fraction of observed successes in each cell as the response, and weighting each summand in the log-likelihood calculation by the number of observations in that cell.

$$\sum_in_i \left(\log(p_i)\sum_jy_{ij}/n_i + \log(1 - p_i)(1 - \sum_j(y_{ij}/n_i)\right)$$

This is exactly equal to the model deviance we calculated above, which you can see by pulling in the sum over $j$ in the long form equation as far as possible.

Meanwhile the saturated deviance is different. Since we no longer have 0-1 responses, even with one parameter per observation we can't get exactly 0. Instead the saturated model log-likelihood is

$$\sum_i n_i\left(\log(\sum_jy_{ij}/n_i)\sum_jy_{ij}/n_i + \log(1 - \sum_jy_{ij}/n_i)(1-\sum_jy_{ij}/n_i)\right).$$

In your example, you can verify that twice this amount is the difference between the reported null and residual deviance values for both models.

ni = dfShort$nReps
yavg = dfShort$mortalityP
sum.terms <-ni*(log(yavg)*yavg + log(1 - yavg)*(1 - yavg))
# Need to handle NaN when yavg is exactly 0
sum.terms[1] <- log(1 - yavg[1])*(1 - yavg[1])

2*sum(sum.terms)
fitShortP$deviance - fitLong$deviance
$\endgroup$
  • $\begingroup$ I think you'll have to clarify the expression for deviance of the saturated models. Log of 0 doesn't operate so well. $\endgroup$ – AdamO Jan 12 '18 at 18:58
  • $\begingroup$ Thanks, I should have clarified what I meant. I added an edit to clarify that by 0log(0) I mean 0 in this context. $\endgroup$ – Jonny Lomond Jan 12 '18 at 19:22
  • $\begingroup$ OK but I'm properly confused (forgive me I never covered deviance in any great detail): if you have log(y)y - log(1-y)(1-y) as the saturated model deviance, isn't every observation just 0? $\endgroup$ – AdamO Jan 12 '18 at 19:28
  • 2
    $\begingroup$ The "saturated model" is an imagined model with one parameter per observation. So its predicted probability for every observation is 1 or 0, depending on the actual observed value. So in this case the log likelihood of the saturated model is indeed 0, the data are the only data which could be generated by the saturated model. $\endgroup$ – Jonny Lomond Jan 12 '18 at 19:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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