cómo replicar el comando "margins, atmeans" de stata en R con la biblioteca margins

Question

Estoy tratando de replicar la salida de margins female, atmeans en R que se muestra aquí:

https://stats.idre.ucla.edu/stata/dae/using-margins-for-predicted-probabilities/

puedo recrear la misma configuración mostrada en la página de la ucla con

library(foreign)
library(margins)

x <- read.dta( "https://stats.idre.ucla.edu/stat/data/hsbdemo.dta" )
x$honors <- as.numeric( x$honors == 'enrolled' )

y este código coincide con el logit honores i.female leer línea en la página de ucla

this_model <- glm( honors ~ female + read , data = x , family = binomial() )
summary(this_model)

a partir de aquí, estoy confundido acerca de cómo modificar

summary(margins(this_model))

para que yo reproduzca el stata atmeans mi objetivo es reproducir la salida mostrada en la página de ucla, pero en R en lugar de stata:

margins female, atmeans

Adjusted predictions                              Number of obs   =        200
Model VCE    : OIM

Expression   : Pr(honors), predict()
at           : 0.female        =        .455 (mean)
               1.female        =        .545 (mean)
               read            =       52.23 (mean)

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      female |
          0  |   .1127311   .0350115     3.22   0.001     .0441097    .1813524
          1  |   .2804526   .0509114     5.51   0.000     .1806681    .3802371
------------------------------------------------------------------------------

el R margins páginas de ayuda de la biblioteca discute el atmeans pero no es obvio para mí cómo implementar esto:

atmeans: calcula los efectos marginales en la media (MEMs) de un conjunto de datos en lugar del comportamiento por defecto de calcular los efectos marginales medios (AMEs)

El argumento atmeans de Stata no está implementado en margins() por varias razones, entre ellas porque es posible conseguir el efecto manualmente mediante una operación como data$var <- mean(data$var, na.rm = TRUE) y pasando el marco de datos modificado a margins(x, data = data) .

Answer 1

La función margin sólo informa del efecto.

> mar <-  margins(this_model, 
+                          data = x,
+                          type = "response",
+                          var = c("female"),
+                          at = list(read = mean(x$read))
+ )
> summary(mar)
       factor    read    AME     SE      z      p  lower  upper
 femalefemale 52.2300 0.1677 0.0585 2.8673 0.0041 0.0531 0.2824

Puedes ver que esto es más o menos equivalente al código de stata, por ejemplo $0.2804526-0.1127311 = 0.1677215 \approx 0.1677$

---

Esto se menciona aquí

https://cran.r-project.org/web/packages/margins/vignettes/Introduction.html

La funcionalidad del comando de Stata para producir márgenes predictivos no se no se ha portado, ya que se puede obtener fácilmente del paquete de predicción

---

Así que puedes utilizar la función de predicción. Por ejemplo:

newdata <- data.frame(female = c("male","female"),
                        read = rep(mean(x$read),2))
p <- predict(object = this_model,
             newdata = newdata,
             type = "response",
             se.fit = TRUE)
mult <- qnorm(0.5*(1-0.95))
out <- cbind(p$fit,
             p$se.fit,
             p$fit+p$se.fit*mult,
             p$fit-p$se.fit*mult)
rownames(out) <- levels(newdata$female)[newdata$female]
colnames(out) <- c("margin", "Std.Err", "lower 95% conf", "upper 95% conf")

con salida

> out
          margin    Std.Err lower 95% conf upper 95% conf
male   0.1127310 0.03501142     0.04410991      0.1813522
female 0.2804526 0.05091129     0.18066832      0.3802369

En lo anterior, los intervalos de confianza se calculan manualmente, ya que el predict no lo hace para los objetos glm. Creo que es un poco difícil a calcular el intervalo de confianza para los objetos glm, ya que la estimación de la varianza (¿y quizás también de los grados de libertad efectivos?) no es sencilla. El método anterior, que coincide con la salida de stata, utiliza un cálculo de intervalos de confianza del 95% basado en $\pm 1.96$ veces el error estándar y se basa en una distribución normal (para la estimación de la media/margen con varianza conocida). Esto funciona porque los datos se modelan con una distribución binomial (donde la estimación de la media y la varianza están relacionadas). Pero no funciona en general para el MLG y a veces hay que utilizar una distribución t que da un intervalo ligeramente más amplio (o mucho más grande cuando el tamaño de la muestra es pequeño).

Como alternativa, existe la glm.predict que calcula todo esto basándose en una simulación (no sé qué tipo de simulación).

> library(glm.predict)
> predicts(this_model,"F;mean",sim.count = 1000000, set.seed = 1)
       mean      lower     upper female  read
1 0.1174438 0.06010156 0.2013822   male 52.23
2 0.2832163 0.19213931 0.3900793 female 52.23

Answer 2

siguiendo con la información para replicar la segunda mitad de esta página de ucla:

library(foreign)

x <- read.dta( "https://stats.idre.ucla.edu/stat/data/hsbdemo.dta" )
x$honors <- as.numeric( x$honors == 'enrolled' )
x$female <- as.numeric( x$female == 'female' )

this_model <- glm( honors ~ female + read , data = x , family = binomial() )

newdata <- expand.grid(female = mean(x$female),
                            read = seq(20,70,10))
    p <- predict(object = this_model,
                 newdata = newdata,
                 type = "response",
                 se.fit = TRUE)
    mult <- qnorm(0.5*(1-0.95))
    out <- cbind(p$fit,
             p$se.fit,
                 p$fit+p$se.fit*mult,
                 p$fit-p$se.fit*mult)
rownames(out) <- seq(20,70,10)
colnames(out) <- c("margin", "Std.Err", "lower 95% conf", "upper 95% conf")

coincide con stata margins, at(read=(20(10)70)) atmeans vsquish post de salida:

   margin      Std.Err lower 95% conf upper 95% conf
20 0.002226388 0.001962525   -0.001620091    0.006072867
30 0.009363869 0.006100713   -0.002593308    0.021321046
40 0.038500176 0.016282790    0.006586493    0.070413859
50 0.145023969 0.031199410    0.083874250    0.206173689
60 0.418114778 0.049886309    0.320339409    0.515890147
70 0.752714028 0.067025506    0.621346451    0.884081606