2016-12-22 103 views
2

我在R中運行OLS迴歸,從中得到一對係數。以下是部分代碼:最小二乘迴歸係數非線性函數的標準誤和置信區間

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 

然後我需要計算:(Attacks + Attackslag)/(1 - TreatmentEffectlag)。問題是我可以使用(coeffs[3] + coeffs[4])/(1 - coeffs[2])在R上做這個,但結果是沒有任何p值或置信區間的固定數字,就像計算器會顯示我一樣。

有誰知道是否有任何函數可以用來計算這個置信區間?


編者注

如果目標量是迴歸係數的線性函數,那麼問題降低到一般的線性假設檢驗,其中精確推斷是可能的。

+0

'bootstrap'它。 – user20650

+0

歡迎來到StackOverflow!請閱讀關於[如何提出一個好問題](http://stackoverflow.com/help/how-to-ask)以及如何給出[可重現的示例]的信息(http://stackoverflow.com/questions/ 5963269 /如何對化妝一個偉大-R-重複性,例如/ 5963610)。這會讓其他人更容易幫助你。 – Jaap

+0

你打算如何使用這些?這將決定最佳的迴應。 – Elin

回答

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 
+2

不錯。我打算用delta方法來做這件事。有關此處發生的情況的解釋,可以參考OP至https://ms.mcmaster.ca/~bolker/emdbook/chap7A.pdf的第5部分(具體爲5.3) –

3

下面是使用增量方法從 '汽車' 包的自包含例如:

# 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