プロビットとロジット回帰のRの堅牢でクラスター化された標準誤差

Question

元々、私はStataで非常に直感的なRのクラスター化された標準誤差を持つプロビット/ロジットモデルを実行したいと考えています。ここで私は答えを見つけたLogistic regression with robust clustered standard errors in R。 したがって、StataとRの両方の結果を堅牢な標準エラーとクラスタ化された標準エラーと比較しようとしました。しかし、私はソフトウェア全体の両方の標準エラーの出力がまったく同じではないことに気付きました。しかし、私がここに示唆した方法を使用すればhttps://diffuseprior.wordpress.com/2012/06/15/standard-robust-and-clustered-standard-errors-computed-in-r/。線形回帰のためにRとStataの両方から正確な出力を得ることができます。したがって、私はRで書いたコードが正しくないと思います。ロジットモデルの代わりにプロビットモデルを実行したい場合は、どのコマンドを使用すればよいでしょうか？これを解決するための優雅な選択肢があるかどうか？ありがとう。

Rコード

## 1. linear regression 
library(rms) 
# model<-lm(Sepal.Length~Sepal.Width+Petal.Length+Petal.Width,iris) 
summary(model) 
fit=ols(Sepal.Length~Sepal.Width+Petal.Length+Petal.Width, x=T, y=T, data=iris) 
fit 
robcov(fit) #robust standard error 
robcov(fit, cluster=iris$Species) #clustered standard error 


## 2. logistic regression 
##demo data generation 
set.seed(1234) 
subj<-rep(1:20,each=4) 
con1<-rep(c(1,0),40) 
con2<-rep(c(1,1,0,0),20) 
effect<-rbinom(80,1,0.34) 
data<-data.frame(subj,con1,con2,effect) 
library(foreign);write.dta(data,'demo_data.dta') 

library(rms) 
fit=lrm(effect ~ con1 + con2, x=T, y=T, data=data) 
fit 
robcov(fit) ##robust standard error 
robcov(fit, cluster=data$subj) ## clustered standard error 

Stataのコード

## 1. linear regression 
webuse iris 
reg seplen sepwid petlen petwid 
reg seplen sepwid petlen petwid,r 
reg seplen sepwid petlen petwid,cluster(iris) 


## 2. logistic regression 

use demo_data,clear 
logit effect con1 con2 
logit effect con1 con2,r 
logit effect con1 con2,cluster(subj) 

Answer 1

私は標準誤差を計算するsandwichパッケージを好みます。 1つの理由は、優れたドキュメントです。使用可能なすべてのデフォルトとオプションを明確に示したvignette("sandwich")と、?sandwichとカスタムbreadとmeatを使用する方法を説明しているcorresponding articleを参照してください。

sandwichを使用して、投稿したオプションの違いを把握することができます。差は自由度補正の可能性が高いでしょう。ここでは、単純な線形回帰の比較：

library(rms) 
library(sandwich) 

fitlm <-lm(Sepal.Length~Sepal.Width+Petal.Length+Petal.Width,iris) 

#Your Blog Post: 
X <- model.matrix(fitlm) 
n <- dim(X)[1]; k <- dim(X)[2]; dfc <- n/(n-k)  
u <- matrix(resid(fitlm)) 
meat1 <- t(X) %*% diag(diag(crossprod(t(u)))) %*% X 
Blog <- sqrt(dfc*diag(solve(crossprod(X)) %*% meat1 %*% solve(crossprod(X)))) 

# rms fits: 
fitols <- ols(Sepal.Length~Sepal.Width+Petal.Length+Petal.Width, x=T, y=T, data=iris) 
Harrell <- sqrt(diag(robcov(fitols, method="huber")$var)) 
Harrell_2 <- sqrt(diag(robcov(fitols, method="efron")$var)) 

# Variations available in sandwich:  
variations <- c("const", "HC0", "HC1", "HC2","HC3", "HC4", "HC4m", "HC5") 
Zeileis <- t(sapply(variations, function(x) sqrt(diag(vcovHC(fitlm, type = x))))) 
rbind(Zeileis, Harrell, Harrell_2, Blog) 

      (Intercept) Sepal.Width Petal.Length Petal.Width 
const  0.2507771 0.06664739 0.05671929 0.1275479 
HC0   0.2228915 0.05965267 0.06134461 0.1421440 
HC1   0.2259241 0.06046431 0.06217926 0.1440781 
HC2   0.2275785 0.06087143 0.06277905 0.1454783 
HC3   0.2324199 0.06212735 0.06426019 0.1489170 
HC4   0.2323253 0.06196108 0.06430852 0.1488708 
HC4m  0.2339698 0.06253635 0.06482791 0.1502751 
HC5   0.2274557 0.06077326 0.06279005 0.1454329 
Harrell  0.2228915 0.05965267 0.06134461 0.1421440 
Harrell_2 0.2324199 0.06212735 0.06426019 0.1489170 
Blog  0.2259241 0.06046431 0.06217926 0.1440781 

1. 結果ブログエントリからはHC1と同等です。ブログエントリがStataの出力と似ている場合、StataはHC1を使用します。
2. Frank Harrelの関数は、HC0に似た結果をもたらします。私が理解する限り、これは最初に提案された解決策であり、vignette(sandwich)または?sandwich::vcovHCに記載されている記事を見ると、他の方法はやや優れた特性を持っています。彼らは自由度の調整が異なります。また、robcov(., method = "efron")への呼び出しはHC3に似ています。

いずれの場合でも、同じ出力が必要な場合はHC1を使用するか、適切に分散 - 共分散行列を調整してください。結局のところ、vignette(sandwich)を見て、異なるバージョン間の違いについては、定数で再スケーリングするだけでHC1からHC0になる必要があることがわかります。これはあまり難しくありません。ところで、HC3またはHC4が典型的には、より良い小さなサンプル特性、および影響力のある観察のもとでのそれらの挙動のために好ましいことに留意されたい。したがって、おそらくデフォルト値をStataに変更したいと思うでしょう。

これらの分散共分散行列は、lmtest::coeftestやcar::linearHypothesisなどの適切な関数に渡すことで使用できます。例えば：クラスタ・堅牢な標準誤差については

library(lmtest) 
coeftest(fitlm, vcov=vcovHC(fitlm, "HC1")) 

t test of coefficients: 

       Estimate Std. Error t value Pr(>|t|)  
(Intercept) 1.855997 0.225924 8.2151 1.038e-13 *** 
Sepal.Width 0.650837 0.060464 10.7640 < 2.2e-16 *** 
Petal.Length 0.709132 0.062179 11.4046 < 2.2e-16 *** 
Petal.Width -0.556483 0.144078 -3.8624 0.0001683 *** 
--- 
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

、あなたはサンドイッチの肉を調整する必要があります（?sandwichを参照）、またはそれを行うfunctionを探します。excruciatingdetailhowtodoitwithappropriatecodesorfunctionsにすでにseveralsourcesexplainingがあります。私がここでホイールを再発明する理由はないので、私はこれをスキップします。

リニアモデルと一般化線形モデルでは、クラスタに堅牢な標準誤差を計算する比較的新しい便利なパッケージもあります。 hereを参照してください。