2
$\begingroup$

Using R, I'd like to test whether multiple parameters in a regression model are equal to specific values (by default, are multiple parameters equal to 0).

For example, in this regression model:

score ~ beta0 + beta1*i1 + beta2*i2 + beta3*age + beta4*i1*age + beta5*i2*age

I want to test H0: (beta2 = 0) and (beta5 = 0)

SAS can do this in PROC REG using the TEST statement.

SAS Program editor contents:

proc reg data=tolerate;
  model score=i1 i2 age i1age i2age;
  test i2=0, i2age=0;    * do assc prof have same reg line as asst? ;
run;

SAS Output window contents:

       Test 1 Results for Dependent Variable score
                                Mean
Source             DF         Square    F Value    Pr > F
Numerator           2        0.15581       0.38    0.6859
Denominator        24        0.40678

Here's code in R that starts the analysis:

# data description: http://statacumen.com/teach/ADA2/ADA2_HW_07_S13.pdf
tolerate <- read.csv("http://statacumen.com/teach/ADA2/ADA2_HW_07_tolerate.csv")
tolerate$rank <- factor(tolerate$rank)
tolerate$rank <- relevel(tolerate$rank, "3")
str(tolerate)

tolerate.manual <- data.frame(score = tolerate$score
                                , i1 = (tolerate$rank==1)
                            , i2 = (tolerate$rank==2)
                                , age = tolerate$age
                            , i1age = tolerate$age * (tolerate$rank==1)
                            , i2age = tolerate$age * (tolerate$rank==2)
                            )
lm.man <- lm(score ~ i1 + i2 + age + i1age + i2age, data = tolerate.manual)
summary(lm.man)

I have been unable to find a solution using library multcomp, contrast, or C(). Ideally, I could do this without creating separate terms in the model, but directly from this lm() statement:

lm.s.a.r.ar <- lm(score ~ age*rank, data = tolerate)

This gets close using a Wald test, but I'm looking for the same F-test SAS uses.

summary(lm.s.a.r.ar)
library(aod) # for wald.test()
coef.test.values <- rep(0, length(coef(lm.s.a.r.ar))) # typically, this will be all 0s
wald.test(b = coef(lm.s.a.r.ar) - coef.test.values
        , Sigma = vcov(lm.s.a.r.ar)
        , Terms = c(4,6))

Wald test:
----------
Chi-squared test:
X2 = 0.77, df = 2, P(> X2) = 0.68

Thanks for considering this question.

$\endgroup$
4
$\begingroup$

Do a full and reduced test by fitting a model with all the parameters, then another model with the reduced model (leave out the terms that you want to test with 0). Then do anova(fit1,fit2) to compute the F test.

If you want to test for values other than 0 then use the offset function in the reduced formula.

$\endgroup$
2
$\begingroup$

What Greg said... with exact code following from OP's example:

lm.base <- update(lm.man, ~ . - i2 - i2age)
# or if you prefer to fully specify the formula rather than using update()...
# lm.base <- lm(score ~ i1 + age + i1age, data = tolerate.manual)
anova(lm.base,lm.man)

P.S.

If you tell wald.test() the denominator df (this is stored in lm under $df.residual), it will give you an F statistic. This example follows again from OP's, but is wrapped in a print statement so we can confirm that our p-value matches what has been given above:

print(wald.test(b = coef(lm.s.a.r.ar) - coef.test.values
    , Sigma = vcov(lm.s.a.r.ar)
    , Terms = c(4,6),df=lm.s.a.r.ar$df.residual),digits=4)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.