24 de septiembre de 2014

Inferencia estadística

En esta parte del curso examinaremos las siguientes herramientas de inferencia estadística

  • Cálculo de intervalo de confianza
  • Pruebas de hipótesis para diferencias de medias
    • Pruebas para medias de muestra única
    • Pruebas para medias de dos muestras independientes
    • Pruebas para medias de dos muestras relacionadas

Cargamos los datos de trabajo

Base de datos para estos ejercicios: Familia y roles de género 2012, a descargar de:

http://iop-data.pucp.edu.pe/busqueda/encuesta/71?

Se sugiere descargar también el cuestionario para utilizarlo como referencia de libro de códigos. Descomprimir y grabar el archivo SPSS en el directorio de trabajo de R

# Importar la base de datos del SPSS a un data frame de R
library(foreign)
genero <- as.data.frame(read.spss("IOP_1212_01_B.sav"))
## re-encoding from UTF-8

Preparar las variables de trabajo

En estos ejemplos trabajaremos con las siguientes variables:

  • Edad ideal para que una mujer se case (P1)
  • Edad ideal para que un hombre se case (P2)
  • Edad en la que el entrevistado se casó o empezó a convivir (P23)
  • Horas dedicadas al trabajo doméstico (P19A)
  • Horas que la pareja le dedica al trabajo doméstico (P28A)

Antes de utilizar estas variables es necesario evaluar si necesitan algún tipo de acondicionamiento (identificar valores perdidos o "raros", recodificar, etc)

Acondicionamiento de las variables

Las siguientes transformaciones nos permitirán acondicionar las variables para utilizarlas en el análisis:

genero$p1r <- genero$P1
genero$p2r <- genero$P2
genero$p1r[genero$P1 > 40] <- NA
genero$p2r[genero$P2 == 99] <- NA

genero$p19ar <- genero$P19A
genero$p19ar[genero$P19A >= 140] <- NA

genero$p28ar <- genero$P28A
genero$p28ar[genero$P28A > 120] <- NA

Intervalos de confianza

Cálculo de intervalos de confianza para una media

Para calcular el intervalo de confianza de la media de una variable que proviene de una muestra simple al azar usamos la siguiente fórmula: \[ \begin{aligned} IC_{1-\alpha}=\bar{x} \pm Z\sigma_{\bar{x}} \end{aligned} \]

Donde:

  • \(1-\alpha\) es el nivel de confianza deseado
  • Z es el valor de la distribución normal estándar que representa el nivel de confianza
  • \(\sigma_{\bar{x}}\) es el error estándar de la media que se obtiene mediante:

\[ \begin{aligned} \sigma_{\bar{x}} = \frac{S_{x}}{\sqrt{n}} \end{aligned} \]

Cálculo de intervalo de confianza an 95% para una media

Entonces para calcular el intervalo de confianza para una media podemos proceder de la siguiente manera:

media <- mean(na.omit(genero$p19ar))
desv <- sd(na.omit(genero$p19ar))
N <- length(na.omit(genero$p19ar))
error.est <- desv/sqrt(N)
error <- 2*error.est # Fijamos Z=2 para indicar un nivel de confianza de 95%
lim.inf <- media-error
lim.sup <- media+error
resultado1 <- data.frame(media, desv, N, error.est, error, lim.inf, lim.sup)
resultado1
##   media  desv    N error.est error lim.inf lim.sup
## 1 22.61 19.95 1196    0.5769 1.154   21.46   23.76

Cálculo del intervalo de confianza al 95% usando la distribucion de t

media <- mean(na.omit(genero$p19ar))
desv <- sd(na.omit(genero$p19ar))
N <- length(na.omit(genero$p19ar))
error.est <- desv/sqrt(N)
error <- qt(0.975, df= N-1) * error.est # Usar el cuantil 0.975 de t
lim.inf <- media-error
lim.sup <- media+error
resultado2 <- data.frame(media, desv, N, error.est, error, lim.inf, lim.sup)
resultado2
##   media  desv    N error.est error lim.inf lim.sup
## 1 22.61 19.95 1196    0.5769 1.132   21.48   23.74
resultado1
##   media  desv    N error.est error lim.inf lim.sup
## 1 22.61 19.95 1196    0.5769 1.154   21.46   23.76

Uso de una función ad-hoc para el cálculo del intervalo de confianza

Generamos la función:

int.conf <- function(x, ic = 95) {
  y <- na.omit(x)
  Media <- mean(y)
  Desv.Est <- sd(y)
  N <- length(y)
  Error.Est <- Desv.Est/sqrt(N)
  ci.y <- 1-(((100-ic)/100)/2)
  Error <- Error.Est * qt(ci.y, df = N-1)
  Lim.Inf <- Media - Error
  Lim.Sup <- Media + Error
  result <- data.frame(Media, Desv.Est, N, Error.Est, Error, Lim.Inf, Lim.Sup)
  return(result)
  }

Usamos la función para calcular un intervalo de confianza al 95% (ic por defecto de la función)

int.conf(genero$p19ar)
##   Media Desv.Est    N Error.Est Error Lim.Inf Lim.Sup
## 1 22.61    19.95 1196    0.5769 1.132   21.48   23.74

En vez de un intervalo de 95% podemos usar otros niveles de confianza:

int.conf(genero$p19ar, ic=99) # Intervalo de 99% de confianza
##   Media Desv.Est    N Error.Est Error Lim.Inf Lim.Sup
## 1 22.61    19.95 1196    0.5769 1.488   21.12    24.1
int.conf(genero$p19ar, ic=90) # Intervalo de 90% de confianza
##   Media Desv.Est    N Error.Est  Error Lim.Inf Lim.Sup
## 1 22.61    19.95 1196    0.5769 0.9496   21.66   23.56

Cálculo de intervalo de confianza para una proporción

También podemos calcular un intervalo de confianza para una proporción a partir de una muestra grande. En este caso, la fórmula es parecida a la de la media, lo que cambia es el cálculo del error estándar de la proporción:

\[ \begin{aligned} IC_{1-\alpha}= p \pm Z\sigma_{p} \end{aligned} \]

Donde:

  • \(1-\alpha\) es el nivel de confianza deseado
  • Z es el valor de la distribución normal estándar que representa el nivel de confianza
  • \(\sigma_{p}\) es el error estándar de la propoción que se obtiene mediante:

\[ \begin{aligned} \sigma_{p} = \sqrt{\frac{pq}{n}} \end{aligned} \]

Donde:

  • p es la proporción de la categoría de análisis en nuestra muestra
  • q es igual a 1 - p

Ejemplo

Vamos a usar la variable P32

library(descr)
freq(genero$P32, plot = FALSE)
## genero$P32 
##                                           Frequency Percent Valid Percent
## Hago mucho más de lo que me corresponde         147  12.219        21.304
## Hago algo más de lo que me corresponde          141  11.721        20.435
## Hago más o menos lo que me corresponde          256  21.280        37.101
## Hago algo menos de lo que me corresponde         65   5.403         9.420
## Hago mucho menos de lo que me corresponde        53   4.406         7.681
## No contesta                                      28   2.328         4.058
## NA's                                            513  42.643              
## Total                                          1203 100.000       100.000

Preparamos los datos para el cálculo del intervalo de confianza:

p32r <- na.omit(genero$P32)
cat <- ifelse(p32r=="Hago más o menos lo que me corresponde", 1, 0)
prop.table(table(cat))
## cat
##     0     1 
## 0.629 0.371
p <- mean(cat)
p
## [1] 0.371

Preparamos el resto de la información

n <- length(cat)
error.est.p <- sqrt((p*(1-p))/n)
error.p <- 2 * error.est.p # Usamos Z = 2 para indicar un nivel de confianza del 95%
lim.inf.p <- p - error.p
lim.sup.p <- p + error.p
result.p <- data.frame(p, n, error.est.p, error.p, lim.inf.p, lim.sup.p)
result.p
##       p   n error.est.p error.p lim.inf.p lim.sup.p
## 1 0.371 690     0.01839 0.03678    0.3342    0.4078

Funcion para cálculo de intervalo de confianza de proporción

Generamos la función

int.conf.p <- function(x, cat, ic = 95){
  vcat <- na.omit(x)
  p <- mean(ifelse(vcat==cat, 1, 0))
  n <- length(vcat)
  error.est.p <- sqrt((p*(1-p))/n)
  beta <- 1-(((100-ic)/100)/2)
  error <- error.est.p * qnorm(beta)
  lim.inf.p <- p - error
  lim.sup.p <- p + error
  result.p <- data.frame(p, n, error.est.p, error, lim.inf.p, lim.sup.p)
  return(result.p)
}

Usamos la función

int.conf.p(genero$P32, "Hago más o menos lo que me corresponde")# IC95%
##       p   n error.est.p   error lim.inf.p lim.sup.p
## 1 0.371 690     0.01839 0.03604     0.335    0.4071
int.conf.p(genero$P32, "Hago más o menos lo que me corresponde", ic=99)# IC99%
##       p   n error.est.p   error lim.inf.p lim.sup.p
## 1 0.371 690     0.01839 0.04737    0.3236    0.4184
int.conf.p(genero$P32, "Hago más o menos lo que me corresponde", ic=90)# IC90%
##       p   n error.est.p   error lim.inf.p lim.sup.p
## 1 0.371 690     0.01839 0.03025    0.3408    0.4013

Probamos con otra variable

freq(genero$P10, plot=FALSE)
## genero$P10 
##             Frequency   Percent
## Si               1070  88.94431
## No                132  10.97257
## No contesta         1   0.08313
## Total            1203 100.00000
int.conf.p(genero$P10, "Si")
##        p    n error.est.p   error lim.inf.p lim.sup.p
## 1 0.8894 1203    0.009041 0.01772    0.8717    0.9072

Gráficar el intervalo de confianza

Podemos crear un gráfico del intervalo de confianza de la media de una variable. En este ejemplo usaremos el gráfico de puntos del paquete ggplot2. Vamos a representar gráficamente el promedio de horas semanales dedicadas a labores domésticas que le dedican los hombres y las mujeres.

Primero vamos a generar un data frame que contenga los estadísticos que necesitamos graficar (la media, los límites inferior y superior del intervalo de confianza). Para ello usaremos la función "summarySE" del paquete "Rmisc" (que debe instalarse en su sistema). Esta función summaySE hace algo parecido a la función que hemos creado en una de las diapositivas anteriores.

library(Rmisc)
## Loading required package: lattice
## Loading required package: plyr
df <- summarySE(genero, measurevar="p19ar", groupvars="SEXO", na.rm=T)
df
##        SEXO   N p19ar    sd     se    ci
## 1 Masculino 587  13.7 13.50 0.5574 1.095
## 2  Femenino 609  31.2 21.35 0.8653 1.699

Generar un gráfico de puntos:

Un gráfico de puntos que muestra las medias de p19ar según sexo

library(ggplot2)
graf.punto1 <- ggplot(df, aes(x=SEXO, y=p19ar)) + 
  geom_point(size = 3) + ylim(0, 40)

graf.punto1

Ahora añadimos las barras que indican el intervalo de confianza, con títulos

graf.punto1 +  geom_errorbar(aes(ymin=p19ar-ci, ymax=p19ar+ci), width = 0.2) + 
  ylab("Horas semanales") + 
  ggtitle("Intervalo de confianza al 95% para la media de las horas semanales\n dedicadas al trabajo doméstico, según sexo")

Gráfico de barras para intervalo de confianza

Podemos graficar el mismo concepto con barras:

graf.barra1 <- ggplot(df, aes(x=SEXO, y=p19ar)) + 
  geom_bar(width=0.5, fill="white", colour="black", stat="identity") + 
  geom_errorbar(aes(ymin=p19ar-ci, ymax=p19ar+ci), width = 0.2) +
  ylab("Horas semanales") + ggtitle("Intervalo de confianza al 95% para la media de las horas semanales\n dedicadas al trabajo doméstico, según sexo")

graf.barra1

Algunas variaciones

Según ámbito y sexo del entrevistado

df2 <- summarySE(genero, measurevar="p19ar", groupvars=c("SEXO", "Ambito"), 
                 na.rm=T)
df2
##        SEXO          Ambito   N p19ar    sd     se    ci
## 1 Masculino     Lima-Callao 218 14.49 14.27 0.9666 1.905
## 2 Masculino Interior Urbano 246 12.68 12.26 0.7819 1.540
## 3 Masculino  Interior Rural 123 14.34 14.43 1.3007 2.575
## 4  Femenino     Lima-Callao 229 30.16 19.83 1.3107 2.583
## 5  Femenino Interior Urbano 264 29.62 19.59 1.2057 2.374
## 6  Femenino  Interior Rural 116 36.83 26.71 2.4795 4.912

El gráfico de puntos

graf.punto2 <- ggplot(df2, aes(x=Ambito, y=p19ar, colour=SEXO)) +
  geom_point(size = 3) + ylim(0, 50) + 
  geom_errorbar(aes(ymin=p19ar-ci, ymax=p19ar+ci), width = 0.2) + 
  ylab("Horas semanales") + 
  ggtitle("Intervalo de confianza al 95% para la media de las horas semanales\n dedicadas al trabajo doméstico, según ámbito, por sexo")

graf.punto2

El gráfico de barras

graf.bar2 <- ggplot(df2, aes(x=Ambito, y=p19ar, fill=SEXO)) +
  geom_bar(position="dodge", stat = "identity") +
  geom_errorbar(aes(ymin=p19ar-ci, ymax=p19ar+ci), width = 0.2, 
  position=position_dodge(.9)) + ylab("Horas semanales") + 
  ggtitle("Intervalo de confianza al 95% para la media de las horas semanales\n dedicadas al trabajo doméstico, según ámbito, por sexo")

graf.bar2