6.2 Logistic Regression

Author

Affiliation

Vladimir Buskin

Catholic University of Eichstätt-Ingolstadt

Abstract

Modelling categorical (binary) response variables.

Preparation

Script

%>% You can find the full R script associated with this unit here.

Consider the data from Buskin’s (2025)¹ corpus-study on subject pronoun realisation:

¹ The input data can be downloaded from this OSF repository: https://osf.io/qgnms.

# Load libraries
library(tidyverse)
library(broom)
library(DescTools)
library(lmtest)

# Load data
data_pro <- read.csv("INPUT_pronouns.csv", sep = ",", header = TRUE)

# Inspect data
str(data_pro)
head(data_pro)

Target variable:
- Type of pronominal Reference (‘overt’, ‘null’)
Explanatory variables:
- Grammatical Person (‘1.p.’, ‘2.p’, ‘3.p’ as well as the dummy pronouns ‘it’ and ‘there’)
- Register (the text category in the International Corpus of English; ‘S1A’ are informal conversations, whereas ‘S1B’ comprises formal class lessons)
- Variety of English (British English ‘GB’, Singapore English ‘SING’ and Hong Kong English ‘HK’), and
- Referentiality (‘referential’ with an identifiable referent or ‘non-referential’ with no/generic reference)

head(data_pro)

  Reference Person Register Variety Referentiality
1     overt      3      S1A      GB    referential
2     overt      3      S1A      GB    referential
3     overt      3      S1A      GB    referential
4     overt      3      S1A      GB    referential
5     overt      3      S1A      GB    referential
6     overt      3      S1A      GB    referential

table(data_pro$Reference)


 null overt 
  174  4664

Descriptive overview

Show the code

data_pro_stats1 <- data_pro %>% 
  group_by(Variety) %>% 
  count(Reference) %>% 
  mutate(pct = n/sum(n) * 100)

Variety	Reference	n	pct
GB	null	29	1.817043
GB	overt	1567	98.182957
HK	null	75	4.507212
HK	overt	1589	95.492789
SING	null	70	4.435995
SING	overt	1508	95.564005

Show the code

data_pro_stats2 <- data_pro %>% 
  group_by(Variety, Register) %>% 
  count(Reference) %>% 
  mutate(pct = n/sum(n) * 100)

Logistic regression

In contrast to linear regression, logistic regression models a discrete response variable \(Y\) with two outcomes². In the present study, \(Y\) is pronominal Reference and has the outcomes Reference = null and Reference = overt, which represent null and overt subjects, respectively. Dichotomous variables of this kind are also often coded as yes/no or 1/0.

² Logistic regression can also be used for \(\geq 3\) classes by breaking down the response variable into a series of dichotomous variables. This is also known as multinomial logistic regression or softmax regression.

Another difference from linear regression is the output of the model:

In linear regression, we obtain a predicted value for the continuous response variable we’re interested in. For instance, if we’re modelling reaction times, the model will return an estimated mean reaction time (given the predictors).
In logistic regression, the model will return a probability. In variational linguistics, this may correspond to the probability that a speaker will use one syntactic variant versus the other.

A core component of logistic regression is the logistic function. The rationale for using it is that the output of the function will always lie between \(0\) and \(1\), and it will always denote a probability.

The simple logistic model

Assuming a binary response variable \(Y\) with the values 1 and 0 and a single predictor \(X\), the conditional probability \(P(Y = 1 \mid X)\) is then equivalent to the inverse logit in Equation 1.

\[ P(Y = 1 \mid X) = \frac{e^{\beta_0 + \beta_1X}}{1 + e^{\beta_0 + \beta_1X}}. \tag{1}\]

With some algebraic manipulation it can be shaped into a form that is definitely more familiar:

\[ \log\left(\frac{P(Y = 1 \mid X)}{1 - P(Y = 1 \mid X)}\right) = \beta_0 + \beta_1X. \tag{2}\]

The logistic model has several characteristic components. The fraction \(\frac{P(Y = 1 \mid X)}{1-P(Y = 1 \mid X)}\) represents the odds, which stand for to the probability of one outcome (e.g., Reference = null) compared to the other (e.g., Reference = overt). Their logarithmic transformation are the log odds (or logits) of one outcome versus the other.

Understanding log odds

When interpreting the output of a logistic model, note that

positive log odds indicate an increase in \(\frac{P(Y = 1 \mid X)}{1-P(Y = 1 \mid X)}\), whereas
negative log odds indicate a decrease in \(\frac{P(Y = 1 \mid X)}{1-P(Y = 1 \mid X)}\).

In more concrete terms: If we are interested in the probability that the form of pronominal reference is null (our \(Y\)) while taking into account the extra-linguistic context (Register; our \(X\)), the model would then have the general form in Equation 3.

\[ \log\left(\frac{P(\text{Reference} = \text{null} \mid \text{Register})}{1- P(\text{Reference} = \text{null} \mid \text{Register})}\right) = \beta_0 + \beta_1\text{Register} \tag{3}\]

Multiple logistic regression

If more than one predictor is included, the above equations can be expanded so as to take into account \(p\) slopes for \(p\) independent variables \(X_1, X_2, ..., X_p\).

\[ P(Y = 1 \mid X_1, ..., X_p) = \frac{e^{\beta_0 + \beta_1X_1 + ... + \beta_pX_p}}{1 + e^{\beta_0 + \beta_1X_1 + ... + \beta_pX_p}}. \tag{4}\]

Consequently, the log odds correspond to the sum of linear predictors \(\beta_1X_1 + \beta_2X_2 + ...+ \beta_pX_p\) (cf. Equation 5).

\[ \log\left(\frac{P(Y = 1 \mid X_1, ..., X_p)}{1 - P(Y = 1 \mid X_1, ..., X_p)}\right) = \beta_0 + \sum_{i=1}^{p} \beta_i X_i \tag{5}\]

Odds ratios

To assess the strength of an effect, it is instructive to examine the odds ratios that correspond to the model coefficients. Odds ratios (OR) are defined as

\[ OR(X_1) = e^{\beta_1}. \]

Understanding odds ratios

Essentially, the OR describes the ratio between two odds with respect to another independent variable. This is illustrated for Reference given Register below:

\[ \text{OR}(\text{Reference} \mid \text{Register}) = \frac{\frac{P(\text{Reference} = \text{null} \mid \text{Register} = \text{S1A})}{P(\text{Reference} = \text{overt} \mid \text{Register} = \text{S1A})}}{\frac{P(\text{Reference} = \text{null} \mid \text{Register} = \text{S1B})}{P(\text{Reference} = \text{overt} \mid \text{Register} = \text{S1B})}} \]

Read as: ‘The ratio between the probability of a null vs. overt subject in S1A and the probability of a null vs. overt object in S1B’.

Maximum Likelihood Estimation

In contrast to continuous data, the estimation of parameters for discrete response variables is less straightforward in that there is no unique solution. Rather than finding a regression line that minimises the distance to all data points, the default approach of logistic models is to find the parameter values that are most likely, given the data. Hence this procedure is also known as Maximum Likelihood Estimation (MLE).

The model first makes an assumption about the probability distribution of the data. For categorical data, the binomial distribution is a common choice. Equation 6 indicates the corresponding probability mass function, which describes the probability \(\pi\) of observing \(y\) successes in \(k\) independent Bernoulli trials. In other words, if we tossed a coin \(n = 10\) times and observed \(y = 5\) heads (i.e., 5 successes), what is the probability \(\pi\) of a success?

\[ f(y; k; \pi) = \binom{k}{y} \pi^y (1 - \pi)^{n-y} \tag{6}\]

Now, let \(\beta\) denote some parameter of interest (e.g., the slope coefficient of a logistic regression model). Given some observed data, the likelihood of this parameter can be described in terms of the likelihood function \(L(\beta)\) in Equation 7. It assumes \(n\) binomial probability mass functions with trials \(k = 1\) and computes their product. Since the binomial coefficient \(\binom{n}{y}\) is a constant term, it is typically dropped. In essence, we’re multiplying successes \(\pi^{y_i}_i\) with failures \((1 - \pi_i)^{1-y_i}\) for each data point.

\[ L(\beta) = \prod_{i=1}^n \pi^{y_i}_i (1 - \pi_i)^{1-y_i} \tag{7}\]

Within the context of our pronoun data, this can be understood more intuitively as the product of the probability that a subject is a null pronoun with the probability that it is not:

\[ L(\beta) = \prod_{i=1}^n \pi^{\text{null}}_i (1 - \pi_i)^{1-\text{null}} \tag{8}\]

Conventionally, this expression is log-transformed in order to convert the product into a sum because, for one, sums are easier to handle computationally. The log-likelihood function \(\ell(\beta)\) in Equation 9 forms the basis for a variety of goodness-of-fit measures used to evaluate logistic regression models.

\[ \ell(\beta) = \sum_{i=1}^n y_i \log(\pi_i) + (1 - y_i) \log(1 - \pi_i) \tag{9}\]

And in more concrete terms:

\[ \ell(\beta) = \sum_{i=1}^n \text{null}_i \log(\pi_i) + (1 - \text{null}_i) \log(1 - \pi_i) \tag{10}\]

In a next step, all \(\pi_i\)s would be substituted by their respectively probabilities from the right half of Equation 4.

The goal is to find the value that maximises \(\ell(\beta)\), i.e., the maximum likelihood estimator \(\hat{\beta}\). Approximate solutions can be attained via iterative optimisation techniques (e.g., Newton-Ralphson or Gradient Descent). Sometimes the algorithm may fail to find an optimal solution, which R may report as a model convergence error. For further technical details, see Wood (2006: 63-66) or Agresti & Kateri (2022: 291-294).

Workflow in R

Research question and hypotheses

How do the intra- and extra-linguistic variables suggested in the literature affect subject pronoun realisation (Definite Null Instantiation) in British English, Singapore English and Hong Kong English?

Given the significance level \(\alpha = 0.05\), the hypotheses are:

\[ \begin{aligned} H_0: & \quad \text{None of the predictor coefficients deviate from 0}.\\ H_1: & \quad \text{At least one predictor coefficient deviates from 0}. \end{aligned} \]

These can be restated mathematically as:

\[ \begin{aligned} H_0: & \quad \beta_1 = \beta_2 = \cdots = \beta_p = 0 \\ H_1: & \quad \text{At least one } \beta_i \neq 0 \text{ for } i \in \{1, 2, \ldots, p\} \end{aligned} \]

Convert to factors and specify reference levels

The next step involves specifying reference levels for all categorical variables. This step is very important because it will directly impact the parameter estimation procedure and, consequently, influence our interpretation of the model output.

The reference level of the response is usually chosen such that it corresponds to the unmarked or most frequent case. Since overt pronouns are much more common in the data, the reference level of the Reference variable will be set to Reference = overt. This way, the model coefficients will directly represent the probability of the null subject variant (i.e., the special case) given certain predictor configurations.
The predictor levels need to be specified as well. Among other things, we are interested in how the Asian Englishes pattern relative to British English. Therefore, we will define British English as the baseline for comparison.

We will use the following specifications:

Variable	Factor Levels	Preferred Reference level
Register	S1A, S1B	S1A
Variety	GB, SING, HK	GB
Person	1, 2, 3, it, there	3
Referentiality	referential, non-referential	referential

# Store "Reference" as factor
data_pro$Reference <- as.factor(data_pro$Reference)

## Specify reference level (the 'unmarked' case)
data_pro$Reference <- relevel(data_pro$Reference, "overt")

## Print levels
levels(data_pro$Reference)

[1] "overt" "null"

Repeat the procedure for the remaining categorical variables.

Code

# Store "Register" as factor
data_pro$Register <- as.factor(data_pro$Register)

## Specify reference level
data_pro$Register <- relevel(data_pro$Register, "S1A")

# Store "Variety" as factor
data_pro$Variety <- as.factor(data_pro$Variety)

## Specify reference level
data_pro$Variety <- relevel(data_pro$Variety, "GB")

# Store "Person" as factor
data_pro$Person <- as.factor(data_pro$Person)

## Specify reference level
data_pro$Person <- relevel(data_pro$Person, "3")

# Store "Referentiality" as factor
data_pro$Referentiality <- as.factor(data_pro$Referentiality)

## Specify reference level
data_pro$Referentiality <- relevel(data_pro$Referentiality, "referential")

Model fitting

# With (glm); available in base R
# Note the additional "family" argument!
Reference.glm <- glm(Reference ~ Register + Variety + Person + Referentiality, data = data_pro, family = "binomial")

# View model statistics
summary(Reference.glm)


Call:
glm(formula = Reference ~ Register + Variety + Person + Referentiality, 
    family = "binomial", data = data_pro)

Coefficients: (1 not defined because of singularities)
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    -3.4402     0.2221 -15.489  < 2e-16 ***
RegisterS1B                     0.1102     0.1603   0.688  0.49173    
VarietyHK                       0.9856     0.2248   4.385 1.16e-05 ***
VarietySING                     0.9662     0.2272   4.252 2.12e-05 ***
Person1                        -0.9155     0.1800  -5.086 3.65e-07 ***
Person2                        -1.6676     0.2691  -6.198 5.72e-10 ***
Personit                        0.8130     0.2960   2.747  0.00602 ** 
Personthere                    -2.6367     1.0078  -2.616  0.00889 ** 
Referentialitynon-referential       NA         NA      NA       NA    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1498.8  on 4837  degrees of freedom
Residual deviance: 1387.5  on 4830  degrees of freedom
AIC: 1403.5

Number of Fisher Scoring iterations: 7

# Get Pseudo-R^2 (requires the DescTools library)
DescTools::PseudoR2(Reference.glm, "Nagelkerke")

Nagelkerke 
0.08538986

# Likelihood ratio test
lmtest::lrtest(Reference.glm)

Likelihood ratio test

Model 1: Reference ~ Register + Variety + Person + Referentiality
Model 2: Reference ~ 1
  #Df  LogLik Df  Chisq Pr(>Chisq)    
1   8 -693.75                         
2   1 -749.42 -7 111.33  < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Further model diagnostics

Post-hoc evaluation of the logistic regression model is just as important as it is for linear regression. However, several assumptions can be relaxed. While independence of observations, a linear relationship between predictors and response as well as uncorrelated predictors remain essential, the stipulations on the distribution and constant variance of the residuals may be disregarded.

Multicollinearity can be inspected via the vif() function from the car package:

vif(Reference.glm)

                  RegisterS1B                     VarietyHK 
                           NA                            NA 
                  VarietySING                       Person1 
                           NA                            NA 
                      Person2                      Personit 
                           NA                            NA 
                  Personthere Referentialitynon-referential 
                           NA                            NA

The fact that the variance inflation factors cannot be estimated properly is highly problematic, suggesting severe multicollinearity. In fact, checking for linear dependencies via alias() reveals Referentiality = non-referential as the culprit:

alias(Reference.glm)

Model :
Reference ~ Register + Variety + Person + Referentiality

Complete :
                              (Intercept) RegisterS1B VarietyHK VarietySING
Referentialitynon-referential 0           0           0         0          
                              Person1 Person2 Personit Personthere
Referentialitynon-referential 0       0       1        1

Refitting the model Referentiality solves the problem:

Reference.glm2 <- glm(Reference ~ Register + Variety + Person, data = data_pro, family = "binomial")

summary(Reference.glm2)


Call:
glm(formula = Reference ~ Register + Variety + Person, family = "binomial", 
    data = data_pro)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -3.4402     0.2221 -15.489  < 2e-16 ***
RegisterS1B   0.1102     0.1603   0.688  0.49173    
VarietyHK     0.9856     0.2248   4.385 1.16e-05 ***
VarietySING   0.9662     0.2272   4.252 2.12e-05 ***
Person1      -0.9155     0.1800  -5.086 3.65e-07 ***
Person2      -1.6676     0.2691  -6.198 5.72e-10 ***
Personit      0.8130     0.2960   2.747  0.00602 ** 
Personthere  -2.6367     1.0078  -2.616  0.00889 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1498.8  on 4837  degrees of freedom
Residual deviance: 1387.5  on 4830  degrees of freedom
AIC: 1403.5

Number of Fisher Scoring iterations: 7

The VIFs look considerably better now:

vif(Reference.glm2)

RegisterS1B   VarietyHK VarietySING     Person1     Person2    Personit 
   1.036261    2.018748    2.023863    1.129088    1.087312    1.098571 
Personthere 
   1.006525

Removing that problematic variable does not affect the fit at all:

anova(Reference.glm, Reference.glm2)

Analysis of Deviance Table

Model 1: Reference ~ Register + Variety + Person + Referentiality
Model 2: Reference ~ Register + Variety + Person
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1      4830     1387.5                     
2      4830     1387.5  0        0

Stepwise variable selection

With the function drop1(), it is possible to successively remove variables from the complex model to ascertain which ones improve the model significantly (i.e., decrease the deviance and AIC scores).

drop1(Reference.glm, test = "Chisq")

Single term deletions

Model:
Reference ~ Register + Variety + Person + Referentiality
               Df Deviance    AIC    LRT  Pr(>Chi)    
<none>              1387.5 1403.5                     
Register        1   1388.0 1402.0  0.471    0.4926    
Variety         2   1413.7 1425.7 26.179 2.067e-06 ***
Person          3   1471.1 1481.1 83.607 < 2.2e-16 ***
Referentiality  0   1387.5 1403.5  0.000              
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The output suggests that Register and Referentiality do not render the model substantially worse upon their removal.

Confidence intervals and odds ratios

# Tidy the model output
tidy_model <- tidy(Reference.glm2, conf.int = TRUE)

# Remove intercept, compute odds ratios and their CIs
(tidy_model <- tidy_model %>% 
  filter(term != "(Intercept)") %>% 
  mutate(
    odds_ratio = exp(estimate),
    odds.conf.low = exp(conf.low),
    odds.conf.high = exp(conf.high)
  )
)

term	estimate	std.error	statistic	p.value	conf.low	conf.high	odds_ratio	odds.conf.low	odds.conf.high
RegisterS1B	0.1102450	0.1603417	0.6875627	0.4917282	-0.2061670	0.4233684	1.1165516	0.8136972	1.5270967
VarietyHK	0.9856195	0.2247932	4.3845615	0.0000116	0.5564413	1.4409284	2.6794714	1.7444534	4.2246161
VarietySING	0.9661514	0.2272470	4.2515473	0.0000212	0.5315430	1.4257107	2.6278115	1.7015557	4.1608138
Person1	-0.9155357	0.1799959	-5.0864250	0.0000004	-1.2738888	-0.5667152	0.4003021	0.2797417	0.5673861
Person2	-1.6676297	0.2690625	-6.1979274	0.0000000	-2.2280126	-1.1666544	0.1886938	0.1077423	0.3114070
Personit	0.8129708	0.2959570	2.7469221	0.0060157	0.2024969	1.3697902	2.2545961	1.2244562	3.9345253
Personthere	-2.6367446	1.0077634	-2.6164322	0.0088854	-5.5095493	-1.1242383	0.0715940	0.0040479	0.3248998

Visualisation

Plot model coefficients:

Code

# Create the coefficient plot
ggplot(tidy_model, aes(x = estimate, y = term)) +
  geom_point() +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "steelblue") +
  theme_minimal() +
  labs(
    x = "Coefficient Estimate (log-odds)",
    y = "Predictor",
    title = "Coefficient Estimates with Confidence Intervals",
    caption = "*Note that the CIs of singificant predictors do not include 0."
  )

Code

# Plot odds ratios
ggplot(tidy_model, aes(x = exp(estimate), y = term)) +
  geom_point() +
  geom_errorbarh(aes(xmin = odds.conf.low, xmax = odds.conf.high), height = 0.2) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "steelblue") +
  theme_minimal() +
  labs(
    x = "Coefficient Estimate (odds ratios)",
    y = "Predictor",
    title = "Odds ratios with Confidence Intervals",
    caption = "*Note that the CIs of singificant predictors do not include 1."
  )

Interpret the model

The full logistic regression model is significantly better than a null model with no predictors (Likelihood Ratio Test: \(\chi^2 = 111.33\), \(df = 7\), \(p < 0.001\)) and has acceptable fit (Nagelkerke’s-\(R^2\) = \(0.09\)). However, heavy multicollinearity has been detected and traced back to Referentiality; removing this feature resolved the issue and did not affect model fit.

The model coefficients of the reduced model indicate that null subjects are significantly more likely in Singapore English compared to British English (Estimate = 0.97, 95% CI [0.53, 1.43], \(p < 0.001\)). This effect is moderate with an \(OR\) of 2.63 (95% CI [1.70, 4.16]), suggesting that the odds of subject omission are approximately 2.6 times higher in the Singaporean variety.

…

Looking under the hood

Some hidden metrics

Beside the coefficients, the model object Reference.glm2 stores a variety of additional metrics, some of which appear to be rather cryptic:

# A tibble: 1 × 8
  null.deviance df.null logLik   AIC   BIC deviance df.residual  nobs
          <dbl>   <int>  <dbl> <dbl> <dbl>    <dbl>       <int> <int>
1         1499.    4837  -694. 1404. 1455.    1388.        4830  4838

logLik is the log-likelihood of the model at the maximum likelihood estimate \(\ell(\hat{\beta})\) (cf. Equation 9). Closer to 0 (less negative) means a better fit; however, this measure isn’t typically reported.

logLik(Reference.glm2)

'log Lik.' -693.7532 (df=8)

deviance (residual deviance) is \(-2\ell(\hat\beta)\) for the fitted model. It’s the logistic-regression analogue of the residual sum of squares in linear regression. Once again, smaller values indicate better fit.

Reference.glm$deviance

[1] 1387.506

-2*logLik(Reference.glm2)

'log Lik.' 1387.506 (df=8)

null.deviance is the deviance of an intercept-only model with no predictors.

Reference.glm$null.deviance

[1] 1498.836

These two deviance measures can be used to obtain Pseudo-\(R^2\) measure to assess the overall fit of the model.

n <- nrow(Reference.glm2$model)
D_model <- Reference.glm2$deviance
D_null <- Reference.glm2$null.deviance

# Get Cox & Snell's R^2; it's not used widely because of its awkward scale
R2_CS  <- 1 - exp(-(D_null - D_model) / n)

# Equivalent to
PseudoR2(Reference.glm2, "CoxSnell")

 CoxSnell 
0.0227487

# Rescale to [0,1]
R2_max <- 1 - exp(-D_null / n) # get maximum possible value Cox & Snell's R^2 can take on
R2_N   <- R2_CS / R2_max # normalise

# Equivalent to
PseudoR2(Reference.glm2, "Nagelkerke")

Nagelkerke 
0.08538986

AIC/BIC: Akaike Information Criterion/Bayesian Information Criterion. AIC (\(-2\ell(\hat{\beta}) + 2 \times (\text{number of parameters})\)) penalises the model for including too many predictors, especially if they are uninformative; this penalty is even stronger in the BIC. Typically, models with lower AIC/BIC are preferred.

AIC(Reference.glm2)

[1] 1403.506

# Number of parameters = number of estimated coefficients (including intercept)
n_params <- length(coef(Reference.glm2))

# Calculate AIC manually
-2*logLik(Reference.glm2) + 2*n_params

'log Lik.' 1403.506 (df=8)

df.null and df.residual: Degrees of freedom for the fitted and null models respectively; they correspond to (\(n - \text{number of parameters}\)).

# Example: residual df
length(Reference.glm2$y) - n_params

[1] 4830

nobs: Mumber of observations \(n\).

length(Reference.glm2$y)

[1] 4838

How do you actually estimate \(\beta\) in logistic regression?

Recall the likelihood function:

\[ \ell(\beta) = \sum_{i=1}^n y_i \log(\pi_i) + (1 - y_i) \log(1 - \pi_i) \] We need to find the parameter value of \(\beta\) for which the likelihood is as high as possible. This is achieved by setting all derivatives to zero:

\[ \frac{\partial\ell(\beta)}{\partial \beta} = \sum_{i=1}^n x_{ij}(y_i - \pi_i) = 0 \quad \text{ for } j = 0, 1, \dots, p \]

These are known as the score equations. As there is no closed-form solution for this problem, we will need to use an algorithm that will provide us with an approximate answer. A popular approach is the Newton-Raphson method, which requires us to compute a matrix of second partial derivatives known as the Hessian:

\[ \frac{\partial^2\ell(\beta)}{\partial\boldsymbol{\beta} \partial\boldsymbol{\beta}^\top} = -\mathbf{X}^\top\mathbf{W}\mathbf{X} \quad \text{where }\mathbf{W} = \text{diag}(\boldsymbol{\pi}_i(1- \boldsymbol{\pi}_i)) \] The algorithm continuously updates its current solution using the following rule until it finds a stable solution:

\[ \beta^{\text{new}} = \beta^{\text{old}} - \left (\frac{\partial^2\ell(\beta)}{\partial\boldsymbol{\beta} \partial\boldsymbol{\beta}^\top}\right)^{-1} \frac{\partial\ell(\beta)}{\partial \beta}, \] usually starting at 0. For illustration, here is a strongly simplified and computationally inefficient implementation in R:

X <- model.matrix(Reference ~ Register + Variety + Person, data = data_pro)
y <- as.numeric(data_pro$Reference == "null")

# Initialise beta at 0
beta <- rep(0, ncol(X))

# Newton-Raphson (also known as "iteratively reweighted least squares")
for (iter in 1:20) {
  
  pi <- plogis(X %*% beta)                     # π_i = p(x_i; beta)
  
  score <- t(X) %*% (y - pi)                   # score eq. X^T(y - p_i)
  
  W <- diag(as.vector(pi * (1 - pi)))           # weights W = diag(π_i(1-π_i))
  
  hessian <- -t(X) %*% W %*% X                  # Hessian -X^T W X
  
  beta <- beta - solve(hessian) %*% score       # Newton step
}

	term	manual_coef	glm_coef
(Intercept)	(Intercept)	-3.440203	-3.440203
RegisterS1B	RegisterS1B	0.110245	0.110245
VarietyHK	VarietyHK	0.985620	0.985620
VarietySING	VarietySING	0.966151	0.966151
Person1	Person1	-0.915536	-0.915536
Person2	Person2	-1.667630	-1.667630
Personit	Personit	0.812971	0.812971
Personthere	Personthere	-2.636748	-2.636745

References

Agresti, Alan, and Maria Kateri. 2022. Foundations of Statistics for Data Scientists: With r and Python. Boca Raton: CRC Press.

Buskin, Vladimir. 2025. “Definite Null Instantiation in English(es): A Usage-based Construction Grammar approach.” Constructions and Frames. https://doi.org/10.1075/cf.24007.bus.

Hosmer, David W., and Stanley Lemeshow. 2008. Applied Logistic Regression. 2nd ed. New York: Wiley.

James, Gareth, Daniela Witten, Trevor Hastie, and Robert Tibshirani. 2021. An Introduction to Statistical Learning: With Applications in r. New York: Springer. https://doi.org/10.1007/978-1-0716-1418-1.

Levshina, Natalia. 2015. How to Do Linguistics with r: Data Exploration and Statistical Analysis. Amsterdam; Philadelphia: John Benjamins Publishing Company.

Winter, Bodo. 2020. Statistics for Linguists: An Introduction Using r. New York; London: Routledge.

Wood, Simon N. 2006. Generalized Additive Models: An Introduction with R. Boca Raton: Chapman & Hall/CRC.

Recommended reading

Preparation

Descriptive overview

Logistic regression

The simple logistic model

Multiple logistic regression

Odds ratios

Maximum Likelihood Estimation

Workflow in R

Research question and hypotheses

Convert to factors and specify reference levels

Model fitting

Further model diagnostics

Confidence intervals and odds ratios

Visualisation

Interpret the model

Looking under the hood

Some hidden metrics

How do you actually estimate \(\beta\) in logistic regression?

References