R linearHypothesis (Car Package) Testing for Superiority and Not Just Equality
I'm Using the Corrosion Dataset from the Faraway Package in R. I Have Built a Simple Linear Model Using: # Load the Data Data(Corrosion) Help(Corrosion) #...
I'm using the corrosion dataset from the faraway package in R.
I have built a simple linear model using:
# Load the data
data(corrosion)
help(corrosion) # Display in the "Help" window some informations
head(corrosion)
# Simple linear regression
model <- lm(loss ~ Fe, data = corrosion)
summary(model)
One of the question is: "You want to test the null hypothesis that the same expected weight loss is equal to 95mg/dm2/day at 1.5% of iron content"
To answer this question I use the linearHypothesis function from the car package.
# H0: C*beta = rhs
linearHypothesis(model, c(1, 1.5), rhs = 95)
And it give me the p-value for this test.
So, now my question is: if null hypothesis H0 is like: "weight loss of at least 95mg/dm2/day", how can I test it?
In other words, the first question the equation of H0 was: Hypothesis: (Intercept) + 1.5 Fe = 95 And now I want to test this H0: Hypothesis: (Intercept) + 1.5 Fe >= 95
Thanks in advance for your help!
2 Answers
The answer based on t_value <- summary(model)$coefficients["Fe", "t value"] is incorrect. That t_value would be used in a test of H0: Fe = 0, not in a test of H0: Hypothesis: (Intercept) + 1.5 Fe >= 95.
The difficulty here is that the linearHypothesis() function doesn't do one-sided tests.
One way to do the one-sided test is to use the relation between p-values of one-sided and two-sided tests that works in linear models when the parameter estimates have normally distributed errors. The relation is as follows:
If the test of H0: param = x versus H1: param != x has p-value p, then the test of H0: param >= x versus H1: param < x will have p-value of either p/2 or 1-p/2. It will be p/2 when the estimate of the parameter is less than x, and 1-p/2 when it is greater.
In your case, the parameter is (intercept) + 1.5*Fe whose estimate can be found from summary(model) to be 129.787 - 1.5*24.02 = 93.757, which is less than 95, so you should use p/2 = 0.3097/2 = 0.15485.
Edited to add: The answer above is really an answer to the statistical question, rather than the programming question. Here's how to do the programming:
library(faraway)
library(car)
#> Loading required package: carData
#>
#> Attaching package: 'car'
#> The following objects are masked from 'package:faraway':
#>
#> logit, vif
# Load the data
data(corrosion)
# Simple linear regression
model <- lm(loss ~ Fe, data = corrosion)
# Two-sided test of H0: (intercept) + 1.5 Fe = 95
test <- linearHypothesis(model, c(1, 1.5), rhs = 95)
test
#> Linear hypothesis test
#>
#> Hypothesis:
#> (Intercept) + 1.5 Fe = 95
#>
#> Model 1: restricted model
#> Model 2: loss ~ Fe
#>
#> Res.Df RSS Df Sum of Sq F Pr(>F)
#> 1 12 113.45
#> 2 11 102.85 1 10.603 1.1341 0.3097
p <- test$"Pr(>F)"[2]
p
#> [1] 0.3097301
estimate <- coef(model) %*% c(1, 1.5)
estimate
#> [,1]
#> [1,] 93.75676
onesided <- if (estimate < 95) p/2 else 1 - p/2
onesided
#> [1] 0.1548651
Created on 2023-03-05 with reprex v2.0.2
Adrien Riaux, One wayout could be...
First, extract t-value from your fitted model:
# Simple linear regression
model <- lm(loss ~ Fe, data = corrosion)
To test the null hypothesis H0: (intercept) + 1.5 Fe >= 95, we need to calculate the t-value using the estimated coefficients and their standard errors:
# first, extract the estimated coefficients and their standard errors
b <- coef(model)
SE_of_b <- sqrt(diag(vcov(model)))
# the, calculate the t-value for `(intercept) + 1.5 Fe = 95`
t_value <- (b[1] + 1.5 * b[2] - 95) / sqrt(SE_of_b[1]^2 + (1.5 * SE_of_b[2])^2)
And then, calculate the p-value:
# finally, calculate the 'p-values' (for the one-sided test where the
# null hypothesis is '(intercept) + 1.5 Fe >= 95')
p_value <- pt(t_value, df = nrow(corrosion) - 2, lower.tail = TRUE)
Finally,
# if the p-value is less than 0.5 (significance level) and t_value is less than 0,
# reject Null Hypothesis (H0), otherwise, we fail to reject H0
ifelse((p_value < 0.5 & t_value < 0), "Reject H0: loss is less than 95mg/dm2/day", "Failed to reject the H0: weight loss is at least 95mg/dm2/day")
Here, I am assuming the null hypothsis to be defined as:
Null hypothesis (H0): (intercept) + 1.5*Fe >= 95
and,
Alternative hypothesis: (intercept) + 1.5*Fe < 95
Let me know if this gives you some ideas...
Sarah Jenkins is a veteran tech journalist with over 12 years of experience covering artificial intelligence, mobile innovations, and digital ethics. Her insights have appeared in leading technology publications worldwide.
