R linearHypothesis (car package) testing for superiority and not just equality

578 Views Asked by At

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

There are 2 best solutions below

2
On BEST ANSWER

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

3
On

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...