Return to Question

For the best-known GLM, the normal linear model, the three types of inference provide identical results.

I tried this out in R to see what would happen, and I got different p-values when I did my own likelihood ratio test versus the default printout in "summary()" that uses Wald, so something about my interpretation of Agresti's comment is incorrect.

Agresti, Alan. Foundations of linear and generalized linear models. John Wiley & Sons, 2015.

For the best-known GLM, the normal linear model, the three types of inference provide identical results.

I tried this out in R to see what would happen, and I got different p-values when I did my own likelihood ratio test versus the default printout in summary() that uses Wald, so something about my interpretation of Agresti's comment is incorrect.

Agresti, Alan. Foundations of Linear and Generalized Linear Models. John Wiley & Sons, 2015.

In Foundations of Linear and Generalized Linear Models, Agresti makes a comment on page 131 about likelihood ratio, Wald, and Score testing of regression parameters.

In Foundations of Linear and Generalized Linear Models, Agresti makes a comment on page 131 about likelihood ratio, Wald, and Score testing of regression parameters.

set.seed(2020)
N <- 100
x <- rbinom(N, 1, 0.5)
err <- rnorm(N)
y <- 0.5*x + err
G0 <- glm(y~1, family="gaussian")
G1 <- glm(y~x, family="gaussian")
test_stat <- summary(G0)$deviance - summary(G1)$deviance
df <- dim(summary(G1)$coefficients)[1] - dim(summary(G0)$coefficients)[1]
p.value <- 1-pchisq(test_stat, df)
p.value
summary(G1)$coefficients[2,4]

set.seed(2020)
N <- 100 # sample size
R <- 1000 # number of simulations
alpha <- 0.05
lrt_r <- wld_r <- rep(0,R)
for (i in 1:R){
    x <- rbinom(N, 1, 0.5)
    err <- rnorm(N)
    y <- 0.5*x + err
    G0 <- glm(y~1, family="gaussian") # intercept-only model
    G1 <- glm(y~x, family="gaussian") # model with x as a predictor
    test_stat <- summary(G0)$deviance - summary(G1)$deviance
    df <- dim(summary(G1)$coefficients)[1] - dim(summary(G0)$coefficients)[1]
    
    lr <- 1-pchisq(test_stat, df) # likelihood ratio test p-value
    wd <- summary(G1)$coefficients[2,4] # Wald test p-value
    
    # check if the p-values warrant rejection at the level of alpha
    #
    if (lr <= alpha){lrt_r[i] <- 1}
    if (wd <= alpha){wld_r[i] <- 1}
}

# Check the power of each test
#
sum(lrt_r)/R*100 # 70.4%
sum(wld_r)/R*100 # 69.9%

set.seed(2020)
N <- 100
x <- rbinom(N, 1, 0.5)
err <- rnorm(N)
y <- 0.5*x + err
G0 <- glm(y ~ 1, family="gaussian")
G1 <- glm(y ~ x, family="gaussian")
test_stat <- summary(G0)$deviance - 
    summary(G1)$deviance
df <- dim(summary(G1)$coefficients)[1] - 
    dim(summary(G0)$coefficients)[1]
p.value <- 1-pchisq(test_stat, df)
p.value
summary(G1)$coefficients[2, 4]

set.seed(2020)
N <- 100 # sample size
R <- 1000 # number of simulations
alpha <- 0.05
lrt_r <- wld_r <- rep(0,R)
for (i in 1:R){
    x <- rbinom(N, 1, 0.5)
    err <- rnorm(N)
    y <- 0.5*x + err
    G0 <- glm(y ~ 1, family="gaussian") 
                # intercept-only model
    G1 <- glm(y ~ x, family="gaussian") 
           # model with x as a predictor
    test_stat <- summary(G0)$deviance - 
    summary(G1)$deviance
    df <- dim(summary(G1)$coefficients)[1] - 
       dim(summary(G0)$coefficients)[1]
    
    lr <- 1-pchisq(test_stat, df) 
        # likelihood ratio test p-value
    wd <- summary(G1)$coefficients[2, 4] 
        # Wald test p-value
    
    # check if the p-values warrant rejection at the level of alpha
    #
    if (lr <= alpha){lrt_r[i] <- 1}
    if (wd <= alpha){wld_r[i] <- 1}
}

# Check the power of each test
#
sum(lrt_r)/R*100 # 70.4%
sum(wld_r)/R*100 # 69.9%