ミネソタ大学のIPUMSデータセットから、1990年の米国国勢調査のデータを分析しようとしていRます。データが重み付けsurveyされているため、パッケージを使用しています。世帯データを取得するだけで（そして、物事を単純にするために個人変数を無視して）、（世帯収入）の平均を計算しようとしています。これを行うために、次のコードを持つ関数を使用して調査デザインオブジェクトを作成しました。hhincomesvydesign()
> require(foreign)
> ipums.household <- read.dta("/path/to/stata_export.dta")
> ipums.household[ipums.household$hhincome==9999999, "hhincome"] <- NA # Fix missing
> ipums.hh.design <- svydesign(id=~1, weights=~hhwt, data=ipums.household)
> svymean(ipums.household$hhincome, ipums.hh.design, na.rm=TRUE)
mean SE
[1,] 37029 17.365
ここまでは順調ですね。ただし、（同じデータセットの異なる部分を対象としたコードをStata使用して）で同じ計算を実行しようとすると、異なる標準エラーが発生します。
use "C:\I\Hate\Backslashes\stata_export.dta"
replace hhincome = . if hhincome == 9999999
(933734 real changes made, 933734 to missing)
mean hhincome [fweight = hhwt] # The code from the link above.
Mean estimation Number of obs = 91746420
--------------------------------------------------------------
| Mean Std. Err. [95% Conf. Interval]
-------------+------------------------------------------------
hhincome | 37028.99 3.542749 37022.05 37035.94
--------------------------------------------------------------
そして、の作者であるこの猫の皮を剥ぐ別の方法を見ると、survey周波数の重み付けについて次のような提案があります。
expanded.data<-as.data.frame(lapply(compressed.data,
function(x) rep(x,compressed.data$weights)))
ただし、このコードを機能させることができないようです。
> hh.dataframe <- data.frame(ipums.household$hhincome, ipums.household$hhwt)
> expanded.hh.dataframe <- as.data.frame(lapply(hh.dataframe, function(x) rep(x, hh.dataframe$hhwt)))
Error in rep(x, hh.dataframe$hhwt) : invalid 'times' argument
私はそれを修正できないようです。これは、この問題に関連している可能性があります。
つまり、要約すると：
Stataとで同じ答えが得られないのはなぜRですか？- どちらが正しいですか（またはどちらの場合も私は何か間違ったことをしていますか）？
- ソリューションが機能していると仮定すると
rep()、それはの結果を複製しStataますか？ - それを行う正しい方法は何ですか？（など）で
plyr実装された関数に限定されるのではなく、答えが任意の計算を行うためにパッケージを使用することを可能にする場合は称賛に値しますsurveysvymean()svyglm()
アップデート
したがって、こことIPUMSから電子メールで受け取った優れたヘルプの後、次のコードを使用して、調査の重み付けを適切に処理しています。将来誰かがこの問題を抱える場合に備えて、ここで説明します。
最初のStataの準備
IPUMSは現在、データをにインポートするためのスクリプトを公開していないため、、、、またはからR開始する必要があります。とりあえず固執します。IPUMSからインポートスクリプトを実行することから始めます。次に、続行する前に、次の変数を追加します。StataSASSPSSStata
generate strata = statefip*100000 + puma
これにより、フォーム240001ごとに一意の整数が作成されますPUMA。最初の2桁は州のfipコード（メリーランド州の場合は24）で、最後の4桁は州PUMAごとに一意のIDです。使用するR場合は、これを実行することも役立つ場合があります
generate statefip_num = statefip * 1
.dtaこれにより、ファイルをインポートしRてラベルを適用し、基になる整数が失われるため、ラベルなしで追加の変数が作成されます。
スタタとsvyset
キースが説明したように、調査のサンプリングはStataを呼び出すことによって処理されsvysetます。
個人レベルの分析には、現在次のものを使用しています。
svyset serial [pweight=perwt], strata(strata)
これにより、重みがに設定されperwt、階層化が上記で作成した変数に設定され、世帯serial番号を使用してクラスタリングが考慮されます。複数年使用している場合は、試してみることをお勧めします
generate double yearserial = year*100000000 + serial
縦方向のクラスタリングも説明します。
世帯レベルの分析（年なし）の場合：
svyset serial [pweight=hhwt], strata(strata)
自明である必要があります（この場合、シリアルは実際には不要だと思いますが）。に置き換えるserialとyearserial、時系列が考慮されます。
でそれを行うR
.dta上で説明した追加の変数を含むファイルをインポートstrataし、個々の文字で分析していると仮定します。
require(foreign)
ipums <- read.dta('/path/to/data.dta')
require(survey)
ipums.design <- svydesign(id=~serial, strata=~strata, data=ipums, weights=perwt)
または世帯レベルで：
ipums.hh.design <- svydesign(id=~serial, strata=~strata, data=ipums, weights=hhwt)
誰かがこれがお役に立てば幸いです。IPUMSのDwin、Keith、Brandonに感謝します。