sturges - tabla de frecuencias en r

Ponderación de frecuencia en R, comparando resultados con Stata (3)

1 y 2) El comentario que citó de Lumley se escribió en 2001 y es anterior a cualquiera de sus trabajos publicados con el paquete de encuestas que solo se ha publicado hace unos años. Probablemente estés usando "pesos" en dos sentidos diferentes. (Lumley describe tres sentidos posibles al principio de su libro). La función de encuesta svydesign utiliza ponderaciones de probabilidad en lugar de ponderaciones de frecuencia. Parece probable que estos no sean realmente pesos de frecuencia sino más bien pesos de probabilidad, dado el tamaño masivo de ese conjunto de datos, y eso significaría que el resultado del paquete de la encuesta es correcto y el resultado de Stata incorrecto. Si no está convencido, entonces el paquete de encuesta ofrece la función as.svrepdesign () con la que el libro de Lumley describe cómo crear un vector de peso replicado a partir de un objeto svydesign.

3) Creo que sí, pero como dijo RMN ... "Estaría mal".

4) Ya que está mal (IMO) no es necesario.

Estoy tratando de analizar los datos del conjunto de datos IPUMS de la Universidad de Minnesota para el censo de 1990 en EE . UU . En R Estoy usando el paquete de la survey porque los datos están weighted . Solo tomando los datos del hogar (e ignorando las variables de la persona para mantener las cosas simples), estoy tratando de calcular la media de hhincome ( ingresos del hogar ). Para hacer esto, creé un objeto de diseño de encuesta utilizando la función svydesign() con el siguiente código:

> 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

Hasta ahora tan bueno. Sin embargo, recibo un error estándar diferente si intento el mismo cálculo en Stata (usando el código destinado a una parte diferente del mismo conjunto de datos ):

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 --------------------------------------------------------------

Y, viendo otra forma de pelar a este gato, el autor de la survey , tiene esta sugerencia para ponderar la frecuencia:

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

Sin embargo, parece que no consigo que este código funcione:

> 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

Que parece que no puedo arreglar. Esto puede estar relacionado con este problema .

Así que en suma:

¿Por qué no obtengo las mismas respuestas en Stata y R ?
¿Cuál es correcto (o estoy haciendo algo mal en ambos casos)?
Suponiendo que la solución rep() funcione, ¿replicaría eso los resultados de Stata ?
¿Cuál es la forma correcta de hacerlo? Felicitaciones si la respuesta me permite usar el paquete plyr para realizar cálculos arbitrarios, en lugar de svymean() a las funciones implementadas en svymean() , svyglm() etc.)

Actualizar

Entonces, después de la excelente ayuda que recibí aquí y de IPUMS por correo electrónico, estoy usando el siguiente código para manejar adecuadamente la ponderación de la encuesta. Aquí describo en caso de que alguien más tenga este problema en el futuro.

Preparación inicial de Stata

Dado que IPUMS actualmente no publica scripts para importar sus datos a R , deberá comenzar desde Stata , SAS o SPSS . Me quedaré con Stata por ahora. Comience ejecutando el script de importación desde IPUMS. Luego antes de continuar agrega la siguiente variable:

generate strata = statefip*100000 + puma

Esto crea un número entero único para cada PUMA de la forma 240001, con los dos primeros dígitos como el código fip del estado (24 en el caso de Maryland) y los últimos cuatro un ID de PUMA que es único por estado. Si vas a usar R , también puede resultarte útil ejecutar esto.

generate statefip_num = statefip * 1

Esto creará una variable adicional sin etiquetas, ya que la importación de archivos .dta en R aplica las etiquetas y pierde los enteros subyacentes.

Stata y `svyset`

Como explicó Keith, Stata maneja el muestreo de la encuesta invocando svyset .

Para un análisis de nivel individual ahora uso:

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

Esto establece la ponderación a perwt , la estratificación a la variable que creamos anteriormente y utiliza el número de serial del hogar para tener en cuenta la agrupación. Si estuviéramos utilizando varios años, podríamos intentar

generate double yearserial = year*100000000 + serial

para tener en cuenta el agrupamiento longitudinal también.

Para análisis a nivel de hogar (sin años):

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

Debería explicarse por sí mismo (aunque creo que en este caso la serie es realmente superflua). Reemplazar la serial por años yearserial tendrá en cuenta una serie de tiempo.

Haciendolo en `R`

Suponiendo que está importando un archivo .dta con la variable de strata adicional explicada anteriormente y analizando en la letra individual:

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

O a nivel del hogar:

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

Espero que alguien encuentre esto útil, y muchas gracias a Dwin, Keith y Brandon de IPUMS.

Adición leve para las personas que no tienen acceso a Stata o SAS; (Pondría esto en comentarios, pero ...) La biblioteca SAScii puede usar el archivo de código SAS para leer los datos descargados de IPUMS. El código a leer en los datos es del doc

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 )

No debes usar pesos de frecuencia en Stata. Eso está bastante claro. Si IPUMS no tiene un diseño de encuesta "complejo", solo puede usar:

mean hhincome [pw = hhwt]

O, por conveniencia:

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

Lo bueno de la segunda opción es que puede usarla para diseños de encuestas más complejos (a través de las opciones en svyset . Luego, puede ejecutar muchos comandos sin tener que escribir [pw ...] todo el tiempo.