rstatistical-testlinearhypothesis

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


Solution

  • 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