Ralph M - 1 year ago 51

R Question

I am interested in calculating estimates and standard errors for linear combinations of coefficients after a linear regression in R. For example, suppose I have the regression and test:

`data(mtcars)`

library(multcomp)

lm1 <- lm(mpg ~ cyl + hp, data = mtcars)

summary(glht(lm1, linfct = 'cyl + hp = 0'))

This will estimate the value of the sum of the coefficients on

`cyl`

`hp`

`lm`

But, suppose I want to cluster my standard errors, on a third variable:

`data(mtcars)`

library(multcomp)

library(lmtest)

library(multiwayvcov)

lm1 <- lm(mpg ~ cyl + hp, data = mtcars)

vcv <- cluster.vcov(lm1, cluster = mtcars$am)

ct1 <- coeftest(lm1,vcov. = vcv)

`ct1`

`am`

`ct1`

`glht`

Error in modelparm.default(model, ...) :

no ‘coef’ method for ‘model’ found!

Any advice on how to do the linear hypothesis with the clustered variance covariance matrix?

Thanks!

Answer Source

`glht(ct1, linfct = 'cyl + hp = 0')`

won't work, because `ct1`

is not a `glht`

object and can not be coerced to such via `as.glht`

. I don't know whether there is a package or an existing function to do this, but this is not a difficult job to work out ourselves. The following small function does it:

```
LinearCombTest <- function (lmObject, vars, .vcov = NULL) {
## if `.vcov` missing, use the one returned by `lm`
if (is.null(.vcov)) .vcov <- vcov(lmObject)
## estimated coefficients
beta <- coef(lmObject)
## sum of `vars`
sumvars <- sum(beta[vars])
## get standard errors for sum of `vars`
se <- sum(.vcov[vars, vars]) ^ 0.5
## perform t-test on `sumvars`
tscore <- sumvars / se
pvalue <- 2 * pt(abs(tscore), lmObject$df.residual, lower.tail = FALSE)
## return a matrix
matrix(c(sumvars, se, tscore, pvalue), nrow = 1L,
dimnames = list(paste0(paste0(vars, collapse = " + "), " = 0"),
c("Estimate", "Std. Error", "t value", "Pr(>|t|)")))
}
```

Let's have a test:

```
data(mtcars)
lm1 <- lm(mpg ~ cyl + hp, data = mtcars)
library(multiwayvcov)
vcv <- cluster.vcov(lm1, cluster = mtcars$am)
```

If we leave `.vcov`

unspecified in `LinearCombTest`

, it is as same as `multcomp::glht`

:

```
LinearCombTest(lm1, c("cyl","hp"))
# Estimate Std. Error t value Pr(>|t|)
#cyl + hp = 0 -2.283815 0.5634632 -4.053175 0.0003462092
library(multcomp)
summary(glht(lm1, linfct = 'cyl + hp = 0'))
#Linear Hypotheses:
# Estimate Std. Error t value Pr(>|t|)
#cyl + hp == 0 -2.2838 0.5635 -4.053 0.000346 ***
```

If we provide a covariance, it does what you want:

```
LinearCombTest(lm1, c("cyl","hp"), vcv)
# Estimate Std. Error t value Pr(>|t|)
#cyl + hp = 0 -2.283815 0.7594086 -3.00736 0.005399071
```

**Remark**

`LinearCombTest`

is upgraded at Get p-value for group mean difference without refitting linear model with a new reference level, where we can test any combination with combination coefficients `alpha`

:

```
alpha[1] * vars[1] + alpha[2] * vars[2] + ... + alpha[k] * vars[k]
```

rather than just the sum

```
vars[1] + vars[2] + ... + vars[k]
```