The statistical model equivalent to this R formula

Question

This is how my model is written in R:

glm(formula = prop ~ A * B * C * D , family = quasibinomial, data = data, 
    weights = w)

This is a quasibinomial generalized linear model used on the logit scale regressing on a saturated, full factorial model.

How do you write the statistical notation of this model?

Answer 1

There is a cool package mentioned in the GLMM FAQs called equatiomatic.

Reproducible example:

remotes::install_github("datalorax/equatiomatic")
#> Skipping install of 'equatiomatic' from a github remote, the SHA1 (39e8c3ea) has not changed since last install.
#>   Use `force = TRUE` to force installation
library(equatiomatic)

set.seed(123)

# Number of observations
n <- 1000

# Generate random binary data with factor predictors
data <- data.frame(
  A = as.factor(sample(c("beets", "peas"), replace = TRUE, size = n)),
  B = as.factor(sample(c("honey", "milk"), replace = TRUE, size = n)),
  C = as.factor(sample(c("oats", "nuts"), replace = TRUE, size = n)),
  D = as.factor(sample(c("apples", "pears"), replace = TRUE, size = n))
)

# Generate binary response variable using rbinom
data$prop <- rbinom(n, size = 1, prob = 0.5)  # Probability of success is 0.5 for simplicity
head(data)
#>       A     B    C      D prop
#> 1 beets honey oats  pears    0
#> 2 beets  milk oats apples    1
#> 3 beets  milk oats  pears    1
#> 4  peas honey nuts  pears    1
#> 5 beets  milk nuts apples    0
#> 6  peas honey oats  pears    0

# Fit model
m <- glm(formula = prop ~ A * B * C * D , family = quasibinomial, data = data)

# Display the model summary
summary(m)
#> 
#> Call:
#> glm(formula = prop ~ A * B * C * D, family = quasibinomial, data = data)
#> 
#> Coefficients:
#>                          Estimate Std. Error t value Pr(>|t|)  
#> (Intercept)               -0.2097     0.2477  -0.847   0.3973  
#> Apeas                      0.1746     0.3643   0.479   0.6317  
#> Bmilk                      0.3639     0.3525   1.032   0.3022  
#> Coats                      0.5744     0.3609   1.592   0.1118  
#> Dpears                     0.6006     0.3680   1.632   0.1030  
#> Apeas:Bmilk               -0.3288     0.5090  -0.646   0.5185  
#> Apeas:Coats               -0.6934     0.5143  -1.348   0.1778  
#> Bmilk:Coats               -0.6595     0.5131  -1.285   0.1990  
#> Apeas:Dpears              -0.2601     0.5194  -0.501   0.6166  
#> Bmilk:Dpears              -1.2115     0.5345  -2.267   0.0236 *
#> Coats:Dpears              -0.7364     0.5064  -1.454   0.1462  
#> Apeas:Bmilk:Coats          0.3587     0.7310   0.491   0.6237  
#> Apeas:Bmilk:Dpears         0.7011     0.7436   0.943   0.3460  
#> Apeas:Coats:Dpears         0.2624     0.7208   0.364   0.7159  
#> Bmilk:Coats:Dpears         0.9320     0.7326   1.272   0.2036  
#> Apeas:Bmilk:Coats:Dpears  -0.1278     1.0357  -0.123   0.9018  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for quasibinomial family taken to be 1.01626)
#> 
#>     Null deviance: 1385.8  on 999  degrees of freedom
#> Residual deviance: 1368.2  on 984  degrees of freedom
#> AIC: NA
#> 
#> Number of Fisher Scoring iterations: 4

# Extract equation
equatiomatic::extract_eq(m)

^{Created on 2024-03-04 with reprex v2.1.0}

And that is the output pasted as a formula:

$$ \log\left[ \frac { E( \operatorname{prop} ) }{ 1 - E( \operatorname{prop} ) } \right] = \alpha + \beta_{1}(\operatorname{A}_{\operatorname{peas}}) + \beta_{2}(\operatorname{B}_{\operatorname{milk}}) + \beta_{3}(\operatorname{C}_{\operatorname{oats}}) + \beta_{4}(\operatorname{D}_{\operatorname{pears}}) +\\ \beta_{5}(\operatorname{A}_{\operatorname{peas}} \times \operatorname{B}_{\operatorname{milk}}) + \beta_{6}(\operatorname{A}_{\operatorname{peas}} \times \operatorname{C}_{\operatorname{oats}}) + \beta_{7}(\operatorname{B}_{\operatorname{milk}} \times \operatorname{C}_{\operatorname{oats}}) + \beta_{8}(\operatorname{A}_{\operatorname{peas}} \times \operatorname{D}_{\operatorname{pears}}) +\\ \beta_{9}(\operatorname{B}_{\operatorname{milk}} \times \operatorname{D}_{\operatorname{pears}}) + \beta_{10}(\operatorname{C}_{\operatorname{oats}} \times \operatorname{D}_{\operatorname{pears}}) + \beta_{11}(\operatorname{A}_{\operatorname{peas}} \times \operatorname{B}_{\operatorname{milk}} \times \operatorname{C}_{\operatorname{oats}}) +\\ \beta_{12}(\operatorname{A}_{\operatorname{peas}} \times \operatorname{B}_{\operatorname{milk}} \times \operatorname{D}_{\operatorname{pears}}) + \beta_{13}(\operatorname{A}_{\operatorname{peas}} \times \operatorname{C}_{\operatorname{oats}} \times \operatorname{D}_{\operatorname{pears}}) + \beta_{14}(\operatorname{B}_{\operatorname{milk}} \times \operatorname{C}_{\operatorname{oats}} \times \operatorname{D}_{\operatorname{pears}}) +\\ \beta_{15}(\operatorname{A}_{\operatorname{peas}} \times \operatorname{B}_{\operatorname{milk}} \times \operatorname{C}_{\operatorname{oats}} \times \operatorname{D}_{\operatorname{pears}}) $$

Now comparing Coefficients in summary(m) with the model formula above you can see how it was put together.

When formally describing a statistical model, it is important to address not only how you modeled the expected value but also how you modeled the variance. This is important for understanding the model's structure and assumptions. See @BenBolker's answer for that.

Answer 2

Depending on the venue, you might want to do this:

1. set up any old linear model (needn't be quasibinomial) with the right formula
2. use equatiomatic::extract_eq() to get the (horrible) LaTeX code for the full expression for the linear predictor. (equatiomatic has been archived on CRAN, but you can use remotes::install_github("datalorax/equatiomatic") to get it from GitHub or remotes::install_version("equatiomatic", "0.3.1") to get the last version from CRAN)
3. edit the LaTeX lightly.

dd <- replicate(5, rnorm(10)) |> as.data.frame() |> setNames(c("prop", LETTERS[1:4]))
m = lm(prop ~ A*B*C*D, dd)
equatiomatic::extract_eq(m)

$$ \begin{split} \eta & = \alpha + \beta_{1}(\operatorname{A}) + \beta_{2}(\operatorname{B}) + \beta_{3}(\operatorname{C}) + \beta_{4}(\operatorname{D}) + \beta_{5}(\operatorname{A} \times \operatorname{B}) + \\& \beta_{6}(\operatorname{A} \times \operatorname{C}) + \beta_{7}(\operatorname{B} \times \operatorname{C}) + \beta_{8}(\operatorname{A} \times \operatorname{D}) + \beta_{9}(\operatorname{B} \times \operatorname{D}) + \\& \beta_{10}(\operatorname{C} \times \operatorname{D}) + \beta_{11}(\operatorname{A} \times \operatorname{B} \times \operatorname{C}) + \beta_{12}(\operatorname{A} \times \operatorname{B} \times \operatorname{D}) + \beta_{13}(\operatorname{A} \times \operatorname{C} \times \operatorname{D}) + \\& \beta_{14}(\operatorname{B} \times \operatorname{C} \times \operatorname{D}) + \beta_{15}(\operatorname{A} \times \operatorname{B} \times \operatorname{C} \times \operatorname{D}) \\ \mu & = (1+\exp(-\eta))^{-1} \\ \textrm{Var}(\mu) & = \phi \mu (1-\mu)/w \\ \end{split} $$

• since this is a quasibinomial model, there is no expression like $y \sim \textrm{Binomial}(...)$. The best you can do is say that it's an iteratively reweighted least-squares fit, with specified link and variance functions.
• people sometimes write the first equation as $\textrm{logit}(\mu) = ...$, which (if your audience understands it) avoids the need to specify the inverse-link function explicitly (equation 2)

Answer 3

glm(formula = prop ~ A * B * C * D , family = quasibinomial, data = data, weights = w)

Assume that $A,B,C,D$ are binomial variables with a base case at 0 and alternate at 1.

$$logit(prop) = log\left(\frac{prop}{1-prop}\right) = \beta_0 + \beta_A A + \beta_B B + \beta_C C + \beta_D D + \beta_{AB} AB + \beta_{AC} AC + \beta_{AD} AD + \beta_{BC} BC + \beta_{BD} BD + \beta_{CD} CD + \beta_{ABC} ABC + \beta_{ABD} ABD + \beta_{ACD} ACD + \beta_{BCD} BCD + \beta_{ABCD} ABCD$$

?glm tells you, in the details,

A specification of the form first:second indicates the set of terms obtained by taking the interactions of all terms in first with all terms in second. The specification first*second indicates the cross of first and second. This is the same as first + second + first:second.

?family says this:

The quasibinomial and quasipoisson families differ from the binomial and poisson families only in that the dispersion parameter is not fixed at one, so they can model over-dispersion. For the binomial case see McCullagh and Nelder (1989, pp. 124–8). Although they show that there is (under some restrictions) a model with variance proportional to mean as in the quasi-binomial model, note that glm does not compute maximum-likelihood estimates in that model. The behaviour of S is closer to the quasi- variants.

For the dispersion parameter, see this discussion [1] or McCullagh and Nelder as the help suggests.

The weights enter into the least squares estimates or maximum likelihood estimates.