Cómo conseguir agruparon los valores de p en las pruebas realizadas en varios conjuntos de datos imputados?

Question

El uso de Amelia en R, he obtenido varios conjuntos de datos imputados. Después de eso, he realizado un repetidas medidas de prueba en el programa SPSS. Ahora, yo quiero a la piscina de los resultados de la prueba. Sé que se puede utilizar Rubin reglas (implementado a través de cualquier múltiplo de imputación paquete en R) a la piscina de los medios y los errores estándar, pero ¿cómo puedo piscina valores de p? Es posible? Hay una función en R para hacerlo? Gracias de antemano.

Answer 1

Sí, es posible y, sí, hay R funciones que lo hacen. En lugar de calcular los valores de p de los repetidos análisis con la mano, se puede usar el paquete Zelig, que también es mencionado en la viñeta de la Amelia-paquete (para una más informativo método de ver mi actualización a continuación). Voy a utilizar un ejemplo de la Amelia-viñeta para demostrar esto:

library("Amelia")
data(freetrade)
amelia.out <- amelia(freetrade, m = 15, ts = "year", cs = "country")

library("Zelig")
zelig.fit <- zelig(tariff ~ pop + gdp.pc + year + polity, data = amelia.out$imputations, model = "ls", cite = FALSE)
summary(zelig.fit)

Esta es la salida correspondiente incluyendo $p$-valores:

  Model: ls
  Number of multiply imputed data sets: 15 

Combined results:

Call:
lm(formula = formula, weights = weights, model = F, data = data)

Coefficients:
                Value Std. Error t-stat  p-value
(Intercept)  3.18e+03   7.22e+02   4.41 6.20e-05
pop          3.13e-08   5.59e-09   5.59 4.21e-08
gdp.pc      -2.11e-03   5.53e-04  -3.81 1.64e-04
year        -1.58e+00   3.63e-01  -4.37 7.11e-05
polity       5.52e-01   3.16e-01   1.75 8.41e-02

For combined results from datasets i to j, use summary(x, subset = i:j).
For separate results, use print(summary(x), subset = i:j).

zelig puede caber un host de modelos distintos de mínimos cuadrados.

Para obtener los intervalos de confianza y los grados de libertad para sus estimaciones puede utilizar mitools:

library("mitools")
imp.data <- imputationList(amelia.out$imputations)
mitools.fit <- MIcombine(with(imp.data, lm(tariff ~ polity + pop + gdp.pc + year)))
mitools.res <- summary(mitools.fit)
mitools.res <- cbind(mitools.res, df = mitools.fit$df)
mitools.res

Esto le dará intervalos de confianza y la proporción de la varianza total que es atribuible a la falta de datos:

              results       se    (lower    upper) missInfo    df
(Intercept)  3.18e+03 7.22e+02  1.73e+03  4.63e+03     57 %  45.9
pop          3.13e-08 5.59e-09  2.03e-08  4.23e-08     19 % 392.1
gdp.pc      -2.11e-03 5.53e-04 -3.20e-03 -1.02e-03     21 % 329.4
year        -1.58e+00 3.63e-01 -2.31e+00 -8.54e-01     57 %  45.9
polity       5.52e-01 3.16e-01 -7.58e-02  1.18e+00     41 %  90.8

Por supuesto, usted puede simplemente combine los resultados interesantes en un solo objeto:

combined.results <- merge(mitools.res, zelig.res$coefficients[, c("t-stat", "p-value")], by = "row.names", all.x = TRUE)

Actualización

Después de algo de práctica, he encontrado una manera más flexible para obtener toda la información necesaria mediante la mice-paquete. Para que esto funcione, necesitará modificar el paquete de la as.mids()-función. Uso Gerko la versión publicada en mi pregunta de seguimiento:

as.mids2 <- function(data2, .imp=1, .id=2){
  ini <- mice(data2[data2[, .imp] == 0, -c(.imp, .id)], m = max(as.numeric(data2[, .imp])), maxit=0)
  names  <- names(ini$imp)
      if (!is.null(.id)){
        rownames(ini$data) <- data2[data2[, .imp] == 0, .id]
      }
      for (i in 1:length(names)){
        for(m in 1:(max(as.numeric(data2[, .imp])))){
          if(!is.null(ini$imp[[i]])){
            indic <- data2[, .imp] == m & is.na(data2[data2[, .imp]==0, names[i]])
            ini$imp[[names[i]]][m] <- data2[indic, names[i]]
      }
    } 
  }
  return(ini)
}

Con esto definido, usted puede ir a analizar los conjuntos de datos imputados:

library("mice")
imp.data <- do.call("rbind", amelia.out$imputations)
imp.data <- rbind(freetrade, imp.data)
imp.data$.imp <- as.numeric(rep(c(0:15), each = nrow(freetrade)))
mice.data <- as.mids2(imp.data, .imp = ncol(imp.data), .id = NULL)

mice.fit <- with(mice.data, lm(tariff ~ polity + pop + gdp.pc + year))
mice.res <- summary(pool(mice.fit, method = "rubin1987"))

Esto le dará todos los resultados que obtiene usando Zelig y mitools y más:

                  est       se     t    df Pr(>|t|)     lo 95     hi 95 nmis   fmi lambda
(Intercept)  3.18e+03 7.22e+02  4.41  45.9 6.20e-05  1.73e+03  4.63e+03   NA 0.571  0.552
pop          3.13e-08 5.59e-09  5.59 392.1 4.21e-08  2.03e-08  4.23e-08    0 0.193  0.189
gdp.pc      -2.11e-03 5.53e-04 -3.81 329.4 1.64e-04 -3.20e-03 -1.02e-03    0 0.211  0.206
year        -1.58e+00 3.63e-01 -4.37  45.9 7.11e-05 -2.31e+00 -8.54e-01    0 0.570  0.552
polity       5.52e-01 3.16e-01  1.75  90.8 8.41e-02 -7.58e-02  1.18e+00    2 0.406  0.393

Nota, usando pool() también se puede calcular el $p$-valores con $df$ ajustado para muestras pequeñas por la omisión de la method-parámetro. Lo que es aún mejor, ahora también se puede calcular el $R^2$ y comparar modelos anidados:

pool.r.squared(mice.fit)

mice.fit2 <- with(mice.data, lm(tariff ~ polity + pop + gdp.pc))
pool.compare(mice.fit, mice.fit2, method = "Wald")$pvalue

Answer 2

Normalmente se puede tomar el valor de p mediante la aplicación de Rubin reglas convencionales parámetros estadísticos como la regresión de pesos. Por lo tanto, a menudo no hay necesidad de poner en común los valores de p directamente. También, el cociente de probabilidad estadística pueden ser agrupadas para comparar los modelos. La agrupación de los procedimientos para otras estadísticas se pueden encontrar en mi libro Flexible Imputación de los Datos Faltantes, capítulo 6.

En los casos donde no se conoce la distribución o el método, hay un inédito procedimiento por Licht y Rubin para cara pruebas. He utilizado este procedimiento para piscina p-valores de la wilcoxon() procedimiento, sino que es general y sencilla para adaptarse a otros usos.

Use el procedimiento siguiente SÓLO si todo lo demás falla, por ahora, sabemos muy poco acerca de sus propiedades estadísticas.

lichtrubin <- function(fit){
    ## pools the p-values of a one-sided test according to the Licht-Rubin method
    ## this method pools p-values in the z-score scale, and then transforms back 
    ## the result to the 0-1 scale
    ## Licht C, Rubin DB (2011) unpublished
    if (!is.mira(fit)) stop("Argument 'fit' is not an object of class 'mira'.")
    fitlist <- fit$analyses
        if (!inherits(fitlist[[1]], "htest")) stop("Object fit$analyses[[1]] is not an object of class 'htest'.")
    m <- length(fitlist)
    p <- rep(NA, length = m)
    for (i in 1:m) p[i] <- fitlist[[i]]$p.value
    z <- qnorm(p)  # transform to z-scale
    num <- mean(z)
    den <- sqrt(1 + var(z))
    pnorm( num / den) # average and transform back
}

Answer 3

Usted puede ser que desee comprobar hacia fuera la información sobre p-valor meta-análisis. Un buen punto de partida:

http://en.wikipedia.org/wiki/Fisher%27s_method