2016-12-22 13 views
2

RでOLS回帰を実行していますが、そこから2つの係数が得られます。コードの一部は次のとおりです。最小二乗回帰係数の非線形関数の標準誤差と信頼区間

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/how-to-make-a-great-r-reproducible-example/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

ニース。私はデルタ法でこれを行うつもりでした。ここで何が起こっているのかについては、OPをセクション5(具体的には5.3)https://ms.mcmaster.ca/~bolker/emdbook/chap7A.pdfに紹介することができます –

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 
関連する問題