CODEWITHSUNDEEP
rregressionlinear-regressionlmcoefficients
Gonzalo M.
Gonzalo M.

Reputation: 21

Standard error and confidence interval for nonlinear function of least squares regression coefficients

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.

Upvotes: 2

Views: 547

Answers (2)

Kevin Wright
Kevin Wright

Reputation: 2687

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

Upvotes: 3

Zheyuan Li
Zheyuan Li

Reputation: 73315

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

Upvotes: 4

Related Questions