6

誰かが私がこれらの相対リスク計算（およびその信頼区間）をRで複製するのを手伝ってもらえますか？

ここでは、Stataで使用される同様の手順について説明します。誰かがRでこれを行う方法を教えてもらえますか（私のデータにはクラスターと階層がありますが、より簡単な例を取り上げました）？relrisk.est関数を試しましたが、非常に複雑な設計を処理するため、調査パッケージを使用したいと思います。また、StataとRの推定値を比較したいと思います。ここで提案されているように、ポアソンを使用しています。

###STATA CODE
use http://www.ats.ucla.edu/stat/stata/faq/eyestudy 
tabulate carrot lenses
*same as R binomial svyglm below
xi: glm lenses carrot, fam(bin) 
*switch reference code
char carrot[omit] 1
xi: glm lenses i.carrot, fam(poisson) link(log) nolog robust eform

###R
library(foreign)
library(survey)
D<-read.dta("http://www.ats.ucla.edu/stat/stata/faq/eyestudy.dta")
table(D$lenses,D$carrot)
D$wgt<-rep(1,nrow(D))
Dd<-svydesign(id=~1,data=D,weights=~wgt)
#change category and eform....?
svyglm(lenses~carrot,Dd,family=binomial)
svyglm(lenses~carrot,Dd,family=quasipoisson(log))
4

1 に答える 1

7

例は単純なデータセットであるため、調査パッケージを使用する必要はまったくありません。また、Rを使用して重回帰を学習するときは、より単純な例から始めて、関連する方法の理解を徐々に深めることをお勧めします。

結局のところ、私の意見では、StataとRの哲学は異なります。Stataは、それが何を意味するのか、またはそれがどのように導き出されるのかを知らなくても、簡単に大量の情報をあなたの顔に投げかけます。一方、Rは同じように（またはそれ以上に）強力で用途が広い可能性がありますが、Stataの「自動化」がなく、必要な結果を得るためにゆっくりと進む必要があります。

Stataコードのより直訳は次のとおりです。

require(foreign)
data <- read.dta("http://www.ats.ucla.edu/stat/stata/faq/eyestudy.dta")
with(data, table(carrot, lenses))
glm.out.1 <- glm(lenses ~ carrot, family="binomial", data=data)
summary(glm.out.1)
logLik(glm.out.1)   # get the log-likelihood for the model, as well
glm.out.2 <- glm(lenses ~ factor(carrot, levels=c(1,0)), family="poisson"(link="log"), data=data)
summary(glm.out.2)
logLik(glm.out.2)
exp(cbind(coefficients(glm.out.2), confint(glm.out.2)))

# deriving robust estimates. vcovHC() is in the sandwich package.
require(sandwich)
robSE <- sqrt(diag(vcovHC(glm.out.2, type="HC")))[2]
rob <- coef(glm.out.2)[2]
rob <- exp(c(rob, rob+qnorm(0.025)*robSE, rob-qnorm(0.025)*robSE))
names(rob) <- c("", "2.5 %", "97.5 %")
rob

(link="log")2番目のglm()呼び出しのinは、の場合のデフォルトのリンク関数であるため、必要ないことに注意してくださいfamily="poisson"。

「ロバストな」見積もりを導き出すために、私はこれを読まなければなりませんでした。これは非常に役に立ちました。サンドイッチパッケージの関数を使用してvcovHC()、によって計算されたものとは異なる分散共分散行列を取得しglm()、それを使用して標準誤差とパラメーター推定値を計算する必要があります。

「堅牢な」見積もりは、小数点以下3桁まで、Stataから取得したものとほぼ同じでした。他のすべての結果は完全に同一でした。コードを実行して、自分の目で確かめてください。

ああ、もう1つ、層化されていない設計でglm（）を使用する方法を見つけたら、survey複雑な設計用に変更されたこの分析関数や他の分析関数の対応物を含むパッケージ全体に方法をマッピングします。また、Thomas Lumleyの著書「ComplexSurveys：A guide to analysisusingR」を読むことをお勧めします。

于 2013-01-13T11:43:25.140 に答える