j'essaie d'analyser les données de L'ensemble de données IPUMS de L'Université du Minnesota pour le recensement des États-Unis de 1990R. Je suis l'aide de la survey colis, car les données sont pondéré. En prenant simplement les données du ménage (et en ignorant les variables de la personne pour garder les choses simples), j'essaie de calculer la moyenne de hhincome (revenu du ménage). Pour ce faire j'ai créé un la conception de l'enquête l'objet à l'aide de la svydesign() fonction avec le code suivant:

> 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

jusqu'ici, tout va bien. Cependant, j'obtiens une erreur type différente si je tente le même calcul dans Stata (en utilisant code destiné à une partie différente du même ensemble de données):

use "C:IHateBackslashesstata_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
--------------------------------------------------------------

et, en regardant une autre façon d'écorcher ce chat, l'auteur de survey, a cette suggestion pour la fréquence pondération:

expanded.data<-as.data.frame(lapply(compressed.data,
               function(x) rep(x,compressed.data$weights)))

Cependant, je n'arrive pas à obtenir ce code:

> 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

je n'arrive pas à résoudre. Cela peut être lié à ce problème.

Donc, en somme:

1. Pourquoi ne puis-je pas obtenir les mêmes réponses Stata et R?
2. qui a raison (ou est-ce que je fais quelque chose de mal dans les deux cas)?
3. en supposant que j'ai le rep() la solution de travail, ne serait-ce répliquer Stata'résultats?
4. What is the droit moyen de le faire? Bravo si la réponse me permet d'utiliser le plyr paquet pour faire des calculs arbitraires, plutôt que d'être limité aux fonctions implémentées dans survey (svymean(),svyglm() etc.)

mise à Jour

donc après l'excellente aide que j'ai reçue ici et de IPUMS par courriel, j'utilise le code suivant gérer correctement la pondération des enquêtes. Je décris ici au cas où quelqu'un d'autre a ce problème à l'avenir.

Préparation Initiale Des Statistiques

depuis IPUMS ne publient pas actuellement de scripts pour importer leurs données dans R, vous aurez besoin pour commencer à partir de Stata,SAS, ou SPSS. Je vais rester avec Stata pour l'instant. Commencez par lancer le script d'importation depuis IPUMS. Alors avant de continuer ajouter la variable suivante:

generate strata = statefip*100000 + puma

Cela crée un entier unique pour chaque PUMA du formulaire 240001, les deux premiers chiffres étant le code d'état fip (24 dans le cas du Maryland) et les quatre derniers a PUMA id qui est unique par l'état. Si vous allez utiliser R vous pouvez également trouver utile pour exécuter ce ainsi

generate statefip_num = statefip * 1

cela créera une variable supplémentaire sans étiquettes, puisque importer .dta fichiers dans R appliquer les étiquettes et perdre les entiers sous-jacents.

Stata et svyset

comme Keith l'a expliqué,Stata en invoquant svyset.

Pour un niveau individuel d'analyse maintenant j'utilise:

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

cela fixe la pondération à perwt, la stratification à la variable que nous avons créée ci-dessus, et utilise le ménage serial nombre à prendre en compte clustering. Si nous utilisions plusieurs années, nous pourrions essayer

generate double yearserial = year*100000000 + serial

compte longitudinale de clustering.

pour l'analyse au niveau des ménages (sans années):

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

devrait être explicite (bien que je pense que dans ce cas-ci série est en fait superflu). En remplaçant serialyearserial tiendra compte d'une série chronologique.

en Faisant R

en supposant que vous importez un .dta le fichier avec l' strata variable expliquée ci-dessus et analysée à la lettre individuelle:

require(foreign)
ipums <- read.dta('/path/to/data.dta')
require(survey)
ipums.design <- svydesign(id=~serial, strata=~strata, data=ipums, weights=perwt)

Ou au niveau du ménage:

ipums.hh.design <- svydesign(id=~serial, strata=~strata, data=ipums, weights=hhwt)

espérons que quelqu'un trouve cela utile, et merci beaucoup à Dwin, Keith et Brandon de IPUMS.