Warnings for binomial stats using frequency format

[bwiernik]

With the binomial() family, the response can be specified in three ways (see ?family):

1. Case/Factor: Rows indicating cases, response a factor indicator success/failure
2. Frequency: Rows indicating trial arms, response numeric between 0 and 1 indicating proportion of success, weights indicating number of trials in the arm.
3. Matrix: Rows indicating trial arms, response a two-column matrix, first column indicating number of successes, second column indicating number of failures.

These three specifications yield equivalent coefficients, standard errors, and log-likelihoods, but other overall model fit statistics (e.g., deviance, AIC, df.residual, etc.) aren't accurate because the sample size used for them is wrong. The correct sample sizes for the three cases would be:

1. nobs(model)
2. sum(weights(model))
3. sum(insight::get_response(model))

Because of this, the results for r2() (and insight::n_obs()) aren't accurate when the second or third specifications are used. r2() errors when the third specification is used because it doesn't expect a multi-column response, but it returns inaccurate results silently when the second or third option is used (r2_somers fails for both—frequency form because the Hmisc function rejects nonbinary values; matrix form because it's not expecting a multi-column response).

Would it be possible, when family(model)$family %in% c("binomial", "quasibinomial") to check whether the response is given in frequency or matrix form, then either return an error/warning or instead compute the statistics using appropriate methods?

Example showing errors with `r2_nagelkerke()`

contr <- rbinom(n = 100, size = 1, prob = .2)
treat <- rbinom(n = 100, size = 1, prob = .8)

dat_factor <- data.frame(
  group = factor(c(rep("c", length(contr)), rep("t", length(treat)))),
  outcome = c(contr, treat)
  )
dat_proportion <- data.frame(
  group = factor(c("c", "t")),
  size = c(length(contr), length(treat)),
  outcome_prop = c(sum(contr) / length(contr), sum(treat) / length(treat))
)
dat_matrix <- data.frame(
  group = factor(c("c", "t")),
  size = c(length(contr), length(treat)),
  success = c(sum(contr), sum(treat)),
  failure = c(length(contr) - sum(contr), length(treat) - sum(treat))
)

m_factor <- 
  glm(outcome ~ group,
      family = binomial(link = "logit"),
      data = dat_factor)
m_proportion <- 
  glm(outcome_prop ~ group,
      weights = size,
      family = binomial(link = "logit"),
      data = dat_proportion)
m_matrix <- 
  glm(cbind(success, failure) ~ group,
      family = binomial(link = "logit"),
      data = dat_matrix)

performance::performance::r2_nagelkerke(m_factor)
performance::performance::r2_nagelkerke(m_proportion)
performance::performance::r2_nagelkerke(m_matrix)

[strengejacke]

Is there a safe way to detect which kind of response / model we have?

[strengejacke]

For instance, here we have a matrix response where nobs() returns the correct value:

data(cbpp, package="lme4")
m_test <- glm(cbind(incidence, size-incidence) ~ period, family=binomial, data=cbpp)
nobs(m_test)
#> [1] 56
nrow(cbpp)
#> [1] 56

^{Created on 2020-11-24 by the reprex package (v0.3.0)}

[bwiernik]

That's not the right number though. That's the number of herd-period groups, but not the number of cattle/cases. The R2 values computed based on n = 56 here aren't going to be correct.

library(plyr)
library(tidyr)
library(performance)

dat_agg <- lme4::cbpp
dat_ind <- lme4::cbpp %>%
    rowwise() %>%
    mutate(sick = list(c(rep(0, size - incidence), rep(1, incidence)))) %>%
    ungroup() %>%
    unnest(sick) %>%
    select(-size, -incidence) -> d

m_agg <- glm(cbind(incidence, size-incidence) ~ period, family=binomial, data= dat_agg)
m_ind <- glm(sick ~ period, family=binomial, data= dat_ind)

# same coefficients
parameters::parameters(m_agg)
parameters::parameters(m_ind)

# different R2
r2_tjur(m_agg)
r2_tjur(m_ind)

r2_somers(m_agg)
r2_somers(m_ind)

r2_nagelkerke(m_agg)
r2_nagelkerke(m_ind)

The functions currently only work when the response variable is either a factor or a 0-1 binary variable, unweighted (indicating that the rows are case-level, not group-level). If the response is non-binary (e.g., numeric values greater than 1 or non-integer values) or if it is a matrix, then the data are group-level and should not be used with the current functions.

[strengejacke]

So your answer to this question:

Is there a safe way to detect which kind of response / model we have?

would be:

If the response is non-binary (e.g., numeric values greater than 1 or non-integer values) or if it is a matrix, then the data are group-level and should not be used with the current functions.

And for matrix, nobs would be sum(insight::get_response(model)) , while for non-binary/non-integer, I could assume sum(weights(model))?

[bwiernik]

Yes that is all correct.

[strengejacke]

Ok, this is the solution for the matrix column, however, I think the "minus" is only required when there's a subtraction in the cbind()?

library(dplyr)
library(tidyr)
library(insight)

dat_agg <- lme4::cbpp
dat_ind <- lme4::cbpp %>%
  rowwise() %>%
  mutate(sick = list(c(rep(0, size - incidence), rep(1, incidence)))) %>%
  ungroup() %>%
  unnest(sick) %>%
  select(-size, -incidence) -> d

m_agg <- glm(cbind(incidence, size-incidence) ~ period, family=binomial, data= dat_agg)
m_ind <- glm(sick ~ period, family=binomial, data= dat_ind)

nobs(m_ind)
#> [1] 842
sum(insight::get_response(m_agg)) - sum(lme4::cbpp$incidence)
#> [1] 842

^{Created on 2020-12-04 by the reprex package (v0.3.0)}

[mattansb]

Yes:

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(tidyr)

dat_ind <- lme4::cbpp %>%
  rowwise() %>%
  mutate(sick = list(c(rep(0, size - incidence), rep(1, incidence)))) %>%
  ungroup() %>%
  select(-size, -incidence) %>% 
  unnest(sick)
  
dat_agg <- lme4::cbpp %>% 
  mutate(no_incidence = size - incidence)
  
m_agg1 <- glm(cbind(incidence, size - incidence) ~ period,
              family = binomial,
              data = dat_agg)

m_agg2 <- glm(cbind(incidence, no_incidence) ~ period,
              family = binomial,
              data = dat_agg)

m_ind <- glm(sick ~ period, family = binomial, data = dat_ind)

nobs(m_ind)
#> [1] 842

head(insight::get_response(m_agg1))
#>   incidence size
#> 1         2   14
#> 2         3   12
#> 3         4    9
#> 4         0    5
#> 5         3   22
#> 6         1   18
sum(insight::get_response(m_agg1))
#> [1] 941

head(insight::get_response(m_agg2))
#>   incidence no_incidence
#> 1         2           12
#> 2         3            9
#> 3         4            5
#> 4         0            5
#> 5         3           19
#> 6         1           17
sum(insight::get_response(m_agg2))
#> [1] 842

^{Created on 2020-12-15 by the reprex package (v0.3.0)}

As for option 2:

Frequency: Rows indicating trial arms, response numeric between 0 and 1 indicating proportion of success, weights indicating number of trials in the arm.

@bwiernik do you have any suggestion how we might differentiate between such a case an a "regular" weighted binomial regression? Perhaps when response is a vector and all 0 < response < 1 and weights * response give integers?

[bwiernik]

I think the get_response() behavior is wrong. The response should be a data frame with the first column being successes and the second column be failures. It seems incorrect that the second column in get_response() is returning the values of size rather than size-incidence considering that size-incidence is what was actually specified in the formula.

Compare to the result of model.response():

model.response(model.frame(m_agg))

[bwiernik]

@mattansb Per the documentation for glm(), positive integer weights should be interpreted as frequencies. I think this would only apply when response is a vector and I suppose you could additionally impose the 0 < response < 1 restriction. Note that it must be the "prior" weights that are extracted, not the "working" weights.

Non-NULL weights can be used to indicate that different observations have different dispersions (with the values in weights being inversely proportional to the dispersions); or equivalently, when the elements of weights are positive integers w_i, that each response y_i is the mean of w_i unit-weight observations. For a binomial GLM prior weights are used to give the number of trials when the response is the proportion of successes: they would rarely be used for a Poisson GLM.

[mattansb]

Thanks @bwiernik, this was very informative.

[strengejacke]

I think the get_response() behavior is wrong.

No, I think actually not. See the note in documentation:

Unlike model.frame(), which may contain transformed variables (e.g. if poly() or scale() was used inside the formula to specify the model), get_data() aims at returning the "original", untransformed data (if possible). Consequently, column names are changed accordingly, i.e. "log(x)" will become "x" etc. for all data columns with transformed values.

[bwiernik]

Got it. Didn't realize that. In that case, extracting the response using model.response(model.frame(x)) would probably be a more straightforward approach here versus insight::get_response()