6
$\begingroup$

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?

$\endgroup$
3
  • 1
    $\begingroup$ It would be unhelpful to you for someone to explain what you already know. Can you edit to clarify what you know and where you’re stuck? $\endgroup$
    – Sycorax
    Mar 5 at 0:55
  • $\begingroup$ @Sycorax I don't know how to write the fulll factorial model given the many possible interaction terms. I am also not sure about the quasibinomial. So I think I don't know much and anything would be helpful $\endgroup$
    – ie86
    Mar 5 at 1:00
  • 5
    $\begingroup$ If you model was prop ~ A*B, would you know how to write the model? $\endgroup$
    – Alex J
    Mar 5 at 1:26

3 Answers 3

20
$\begingroup$

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.

$\endgroup$
2
  • 2
    $\begingroup$ This answer is better than mine (especially because you set up the predictors as factors). If you add/adapt the variance function from my answer, I'll delete mine. $\endgroup$
    – Ben Bolker
    Mar 5 at 2:48
  • 2
    $\begingroup$ @BenBolker Not at all! I think our answers complement each other. You also provide other very useful information (in bullet points). And equatiomatic doesn't provide the variance function. So if that's OK with you, I'd like to simply add a note to my answer pointing to the variance function in your post. $\endgroup$
    – Stefan
    Mar 5 at 3:01
14
$\begingroup$

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)
$\endgroup$
1
  • 1
    $\begingroup$ I'm quite sad to hear that equatiomatic is now archived. Even if the output wasn't amazing, it was helpful for quickly extracting the base formula if you had a really complicated model. $\endgroup$ Mar 6 at 3:31
4
$\begingroup$
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.

$\endgroup$

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.