Я пытаюсь проанализировать данные из набора данных IPUMS Университета Миннесоты для переписи населения США 1990 года в R . Я с использованием пакета обследования , потому что данные взвешены . Просто взяв данные о домохозяйстве (и игнорируя личные переменные для простоты), Я пытаюсь вычислить среднее значение хинкома ( семейного дохода ). Для этого я создал объект проектирования , используя функцию 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
--------------------------------------------------------------

И, глядя на другой способ снять шкуру с этой кошки, автор обзора , предлагает следующее предложение для частотного взвешивания:

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

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

Итак, в итоге:

1. Почему не t Я получаю одинаковые ответы в Stata и R ?
2. Какой из них правильный (или я делаю что-то не так в обоих случаях)?
3. Предполагая, что я получил [ Решение 11115333] rep () работает, будет ли оно повторять результаты Stata ?
4. Какой правильный способ сделать это? Престижность, если ответ позволяет мне использовать пакет plyr для выполнения произвольных вычислений, а не ограничиваться функциями, реализованными в обзоре ( svymean () , svyglm () и т. Д.)

Обновление

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

Начальная подготовка Stata

Поскольку 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 , стратификацию для переменной, которую мы создали выше, и использует номер домашнего хозяйства серийный для учета кластеризации. Если бы мы использовали несколько лет, мы могли бы попытаться

generate double yearserial = year*100000000 + serial

также учесть продольную кластеризацию.

Для анализа на уровне домохозяйства (без лет):

svyset serial [pweight=hhwt], strata(strata)

Должно быть понятно (хотя я думаю, что в этом случае серийные номера на самом деле излишни). Замена серийного номера на 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.