Частота взвешивания в R, сравнивая результаты с Stata

Я пытаюсь проанализировать данные из набора данных IPUMS Университета Миннесоты для 1990 США переписи в R. Я использую пакет survey, потому что данные Взвешенный. Просто беря данные домохозяйства (и игнорируя переменные человека, чтобы все было в порядке), я пытаюсь вычислить среднее значение hhincome (доходы домашних хозяйств). Для этого я создал объект дизайна съемки, используя svydesign() со следующим кодом:

> 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

Что я не могу исправить. Это может быть связано с this issue.

Итак, в сумме:

• Почему бы мне не получить одинаковые ответы в Stata и R?
• Какой из них прав (или я делаю что-то неправильно в обоих случаях)?
• Предполагая, что я получил решение rep(), будет ли репликация результатов Stata?
• Что такое правильный способ сделать это? Кудо, если ответ позволяет мне использовать plyr пакет для выполнения произвольных вычислений, а не ограничиваться функции, реализованные в survey (svymean(), svyglm() и т.д.)

Update

Итак, после отличной помощи, которую я получил здесь и из IPUMS по электронной почте, я использую следующий код для правильной обработки взвешивания опроса. Я опишу здесь, если у кого-то еще есть эта проблема в будущем.

Подготовка исходных стат

Поскольку IPUMS в настоящее время не публикуют сценарии для импорта своих данных в R, вам нужно начать с Stata, SAS, или SPSS. На данный момент я придерживаюсь Stata. Начните с запуска импорта script из IPUMS. Затем перед продолжением добавьте следующую переменную:

generate strata = statefip*100000 + puma

Это создает уникальное целое число для каждого PUMA формы 240001 с первыми двумя цифрами в качестве кода состояния fipt (24 в случае Maryland) и последних четырех a PUMA id, который является уникальным для каждого государственной основе. Если вы собираетесь использовать R, вам также может быть полезно также запустить его

generate statefip_num = statefip * 1

Это создаст дополнительную переменную без меток, так как импорт .dta файлов в R применит метки и потеряет основные целые числа.

Stata и 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)

Надеюсь, что кто-то найдет это полезным, и так благодарит Dwin, Keith и Brandon от IPUMS.

Ответы

mean hhincome [pw = hhwt]

svyset [pw = hhwt]
svy: mean hhincome
svy: regress hhincome `x'

library(SAScii)
IPUMS.file.location <- "..\\usa_00007dat\\usa_00007.dat"
IPUMS.SAS.read.in.instructions <- "..\\usa_00007dat\\usa_00007.sas"

#store the IPUMS extract as an R data frame!
IPUMS.df <- 
  read.SAScii ( 
    IPUMS.file.location , 
    IPUMS.SAS.read.in.instructions , 
    zipped = F )