Я пытаюсь проанализировать данные из набора данных 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
Что я не могу исправить.Это может быть связано с этой проблемой .
Итак, в сумме:
- Почему я не получаю одинаковые ответы в
Stata и R? - Какой из них правильный (или я делаю что-то не так в обоих случаях)?
- Предполагается, что у меня работает решение
rep()Будет ли это повторять результаты Stata? - Каков правильный способ сделать это?Престижность, если ответ позволяет мне использовать пакет
plyr для выполнения произвольных вычислений, а не ограничиваться функциями, реализованными в survey (svymean(), svyglm() и т. Д.)
Обновление
Итак, после того, как я получил отличную помощь от IPUMS по электронной почте, я использую следующий код для правильной обработки веса опроса.Я опишу здесь случай, если у кого-то еще возникнет эта проблема в будущем.
Начальная подготовка данных
Поскольку IPUMS в настоящее время не публикует сценарии для импорта своих данных в Rвам нужно будет начать с Stata, SAS или SPSS.Сейчас я буду придерживаться Stata.Начните с запуска сценария импорта из IPUMS.Затем, прежде чем продолжить, добавьте следующую переменную:
generate strata = statefip*100000 + puma
Это создает уникальное целое число для каждого PUMA формы 240001, с первыми двумя цифрами в качестве кода FIP состояния (24 в случае Мэриленда) ипоследние четыре идентификатора PUMA, которые являются уникальными для каждого штата.Если вы собираетесь использовать 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)
Надеюсь, кто-то найдет это полезным, и большое спасибо Двину, Кита и Брэндонамот IPUMS.