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