2

I am running an OLS regression in R from which I get a couple of coefficients. Here's part of the code:

Attacks <- Treat.Terr.Dataset$Attacks[2:30]
Attackslag <- Treat.Terr.Dataset$Attacks[1:29]
TreatmentEffect <- Treat.Terr.Dataset$TreatmentEffect[2:30]
TreatmentEffectlag <- Treat.Terr.Dataset$TreatmentEffect[1:29]

olsreg <- lm(TreatmentEffect ~ TreatmentEffectlag + Attacks + Attackslag)
coeffs<-olsreg$coefficients

And then I would need to calculate: (Attacks + Attackslag) / (1 - TreatmentEffectlag). The problem is that I can do this on R by using (coeffs[3] + coeffs[4]) / (1 - coeffs[2]) but the result is a fixed number without any p-value or confidence interval, just as a calculator would show me.

Does anyone know if there is any function I could use to calculate this with a confidence interval?

Editor note

If the target quantity is a linear function of regression coefficients, then the problem reduces to general linear hypothesis testing where exact inference is possible.

Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248
Gonzalo M.
  • 21
  • 1
  • Welcome to StackOverflow! Please read the info about [how to ask a good question](http://stackoverflow.com/help/how-to-ask) and how to give a [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example/5963610). This will make it much easier for others to help you. – Jaap Dec 22 '16 at 16:40
  • What are you planning to use these for? That will determine the best response. – Elin Dec 22 '16 at 16:50

2 Answers2

4
## variance-covariance of relevant coefficients
V <- vcov(olsreg)[2:4, 2:4]
## point estimate (mean) of relevant coefficients
mu <- coef(olsreg)[2:4]

## From theory of OLS, coefficients are normally distributed: `N(mu, V)`
## We now draw 2000 samples from this multivariate distribution
beta <- MASS::mvrnorm(n = 2000, mu, V)

## With those 2000 samples, you can get 2000 samples for your target quantity
z <- (beta[, 2] + beta[, 3]) / (1 - beta[, 1])

## You can get Monte Carlo standard error, and Monte Carlo Confidence Interval
mean(z)
sd(z)
quantile(z, prob = c(0.025, 0.975))

## You can of course increase sample size from 2000 to 5000
Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248
  • 2
    Nice. I was going to do this by the delta method. For an explanation of what's going on here, you could refer OP to section 5 (specifically 5.3) of https://ms.mcmaster.ca/~bolker/emdbook/chap7A.pdf – Ben Bolker Dec 22 '16 at 16:51
3

Here's a self-contained example using the delta method from the 'car' package:

# Simulate data
dat <- data.frame(Attacks = rnorm(30), Trt=rnorm(30))
dat <- transform(dat, AttacksLag = lag(Attacks), TrtLag = lag(Trt))
dat <- dat[2:30,]

# Fit linear model
m1 <- lm(Trt ~  TrtLag + Attacks + AttacksLag, data=dat)

# Use delta method
require("car")
del1 <- deltaMethod(m1, "(Attacks + AttacksLag) / (1 - TrtLag)")

# Simple Wald-type conf int
del1$Est +  c(-1,1) * del1$SE * qt(1-.1/2, nrow(dat)-length(coef(m1)))
# [1] -0.2921529  0.6723991
Kevin Wright
  • 2,397
  • 22
  • 29