# Práctica 7 de MEE.

### Prevención de riesgos laborales en trabajadores hospitalarios

Francisco García García (2015-06-02)

### A. Introducción

### B. Lectura de los datos

### C. Descripción de variables

### D. Inferencia de proporciones

### E. Inferencia de medias

### F. Inferencia para tablas de contingencia

### G. Bootstrap

## A. Introducción

Se ha realizado un estudio para evaluar el nivel de estrés en el grupo de trabajadores hospitalarios de la provincia de Valencia. Para la obtención de la información, se utilizó un cuestionario para valorar las diferentes exposiciones a riesgos en el trabajo (carga física, carga psíquica y ruido). Además se cuantificó el estrés de cada trabajador.

En el fichero riesgos.txt se presentan los datos de una muestra de esta población de trabajadores. A continuación se describen las variables incluidas:

- id: identificador del trabajador
- edad: edad del trabajador en años
- peso: peso del trabajador en kg.
- genero: hombre, mujer
- sueldo: salario bruto anual expresado en euros
- contrato: tipo de contrato (temporal, fijo)
- jornada: tipo de jornada (completa, parcial)
- carfisica: exposición a carga física, escala 0 (sin exposición) a 5 (mucha exposición)
- carpsiquica: exposición a carga psíquica, escala 0 (sin exposición) a 5 (mucha exposición)
- ruido: exposición ruido, escala 0 (sin exposición) a 5 (mucha exposición)
- estres: nivel de estrés del trabajador, escala 0 (sin estrés) a 10 (mucho estrés)

## B. Lectura de los datos

1) Tras leer la información del archivo riesgo.txt, vamos a conocer cómo están organizados los datos de este grupo de trabajadores:

- ¿Qué nombre reciben las variables de este conjunto de datos?
- ¿Qué dimensión tiene el grupo de datos? (número de columnas y filas, son respectivamente las variables y los trabajadores)
- Veamos el aspecto de los datos para las primeras filas. ¿Y para ver sólo las últimas?
- ¿Cómo visualizamos los datos como si estuviéramos en una hoja de cálculo?

    riesgos <- read.table("riesgos.txt", sep ="\t", as.is = T, header = T)
    colnames(riesgos)

    #>  [1] "id"          "edad"        "peso"        "genero"      "sueldo"
    #>  [6] "contrato"    "jornada"     "carfisica"   "carpsiquica" "ruido"
    #> [11] "estres"

    dim(riesgos)

    #> [1] 91 11

    head(riesgos); tail(riesgos)

    #>   id edad peso genero sueldo contrato  jornada carfisica carpsiquica ruido
    #> 1  1   27   79 hombre  18000     fijo completa         5           3     2
    #> 2  2   33   75 hombre  22000 temporal  parcial         3           3     4
    #> 3  3   37   82 hombre  23000 temporal  parcial         0           4     5
    #> 4  4   32   73 hombre  22000 temporal  parcial         5           4     3
    #> 5  5   31   62 hombre  24000     fijo completa         3           2     2
    #> 6  6   35   77 hombre  23000 temporal  parcial         3           3     2
    #>   estres
    #> 1     10
    #> 2     10
    #> 3     10
    #> 4      2
    #> 5      2
    #> 6      9

    #>    id edad peso genero sueldo contrato  jornada carfisica carpsiquica
    #> 86 86   28   66  mujer  23000     fijo completa         3           2
    #> 87 87   28   67  mujer  23000     fijo completa         1           3
    #> 88 88   34   74  mujer  21000 temporal  parcial         1           5
    #> 89 89   27   59  mujer  23000     fijo completa         3           2
    #> 90 90   28   70  mujer  23000     fijo completa         3           2
    #> 91 91   31   74  mujer  22000 temporal  parcial         3           3
    #>    ruido estres
    #> 86     2      7
    #> 87     0      2
    #> 88     3      2
    #> 89     2      2
    #> 90     1      2
    #> 91     2     NA

    View(riesgos)

## C. Descripción de los datos

2) Descripción de las variables.

Utilizando la función de R: summary, realizamos una rápida descripción de las variables incluidas en la base de datos. Indica cómo son cada una de las variables: categóricas o cuantitativas.

    summary(riesgos)

    #>        id            edad            peso          genero
    #>  Min.   : 1.0   Min.   :23.00   Min.   :58.00   Length:91
    #>  1st Qu.:23.5   1st Qu.:29.00   1st Qu.:67.00   Class :character
    #>  Median :46.0   Median :31.00   Median :72.00   Mode  :character
    #>  Mean   :46.0   Mean   :31.45   Mean   :71.95
    #>  3rd Qu.:68.5   3rd Qu.:34.00   3rd Qu.:76.00
    #>  Max.   :91.0   Max.   :39.00   Max.   :89.00
    #>
    #>      sueldo        contrato           jornada            carfisica
    #>  Min.   :18000   Length:91          Length:91          Min.   :0.0
    #>  1st Qu.:20000   Class :character   Class :character   1st Qu.:2.0
    #>  Median :22000   Mode  :character   Mode  :character   Median :3.0
    #>  Mean   :22165                                         Mean   :3.0
    #>  3rd Qu.:24000                                         3rd Qu.:4.5
    #>  Max.   :26000                                         Max.   :5.0
    #>
    #>   carpsiquica        ruido           estres
    #>  Min.   :2.000   Min.   :0.000   Min.   : 1.000
    #>  1st Qu.:3.000   1st Qu.:2.000   1st Qu.: 2.000
    #>  Median :3.000   Median :3.000   Median : 3.000
    #>  Mean   :3.407   Mean   :2.824   Mean   : 4.756
    #>  3rd Qu.:4.000   3rd Qu.:4.000   3rd Qu.: 8.000
    #>  Max.   :5.000   Max.   :6.000   Max.   :10.000
    #>                                  NA's   :1

3) ¿Qué tipo de distribución presenta las variables estres (nivel de estres de cada trabajador) y edad (edad de cada trabajador)? ¿Podemos asumir que estos datos siguen una distribución normal?

    attach(riesgos)

    boxplot(estres)

    hist(estres)

    qqnorm(estres)
    qqline(estres)

    summary(estres)

    #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's
    #>   1.000   2.000   3.000   4.756   8.000  10.000       1

    boxplot(edad)

    hist(edad)

    qqnorm(edad)
    qqline(edad)

    summary(edad)

    #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    #>   23.00   29.00   31.00   31.45   34.00   39.00

4) ¿Cuántos trabajadores tienen menos de 30 años?

    summary(riesgos)

    #>        id            edad            peso          genero
    #>  Min.   : 1.0   Min.   :23.00   Min.   :58.00   Length:91
    #>  1st Qu.:23.5   1st Qu.:29.00   1st Qu.:67.00   Class :character
    #>  Median :46.0   Median :31.00   Median :72.00   Mode  :character
    #>  Mean   :46.0   Mean   :31.45   Mean   :71.95
    #>  3rd Qu.:68.5   3rd Qu.:34.00   3rd Qu.:76.00
    #>  Max.   :91.0   Max.   :39.00   Max.   :89.00
    #>
    #>      sueldo        contrato           jornada            carfisica
    #>  Min.   :18000   Length:91          Length:91          Min.   :0.0
    #>  1st Qu.:20000   Class :character   Class :character   1st Qu.:2.0
    #>  Median :22000   Mode  :character   Mode  :character   Median :3.0
    #>  Mean   :22165                                         Mean   :3.0
    #>  3rd Qu.:24000                                         3rd Qu.:4.5
    #>  Max.   :26000                                         Max.   :5.0
    #>
    #>   carpsiquica        ruido           estres
    #>  Min.   :2.000   Min.   :0.000   Min.   : 1.000
    #>  1st Qu.:3.000   1st Qu.:2.000   1st Qu.: 2.000
    #>  Median :3.000   Median :3.000   Median : 3.000
    #>  Mean   :3.407   Mean   :2.824   Mean   : 4.756
    #>  3rd Qu.:4.000   3rd Qu.:4.000   3rd Qu.: 8.000
    #>  Max.   :5.000   Max.   :6.000   Max.   :10.000
    #>                                  NA's   :1

    table(riesgos$edad < 30 )

    #>
    #> FALSE  TRUE
    #>    66    25

5) ¿Qué relación existe entre la edad de los trabajadores y su salario?

- Describe gráficamente esta relación y cuantifícala mediante el coeficiente de correlación de Pearson (2 decimales).
- A partir de los resultados obtenidos, ¿qué se concluye?

    plot(riesgos$edad ~ riesgos$sueldo)

    cor.test(riesgos$edad, riesgos$sueldo)

    #>
    #>  Pearson's product-moment correlation
    #>
    #> data:  riesgos$edad and riesgos$sueldo
    #> t = 7.7317, df = 89, p-value = 1.537e-11
    #> alternative hypothesis: true correlation is not equal to 0
    #> 95 percent confidence interval:
    #>  0.4921847 0.7428484
    #> sample estimates:
    #>       cor
    #> 0.6338779

    round(cor.test(riesgos$edad, riesgos$sueldo)$estimate,2)

    #>  cor
    #> 0.63

6) ¿Qué relación existe entre género y jornada? Describe gráfica y numéricamente esta relación mediante:

- Una tabla de contingencia que cuantifique esta relación.
- También un gráfico que muestre como se comportan ambas variables.
- ¿La distribución de la jornada es similar entre hombres y mujeres?

    tabla1 <- table(riesgos$jornada, riesgos$genero)
    tabla1

    #>
    #>            hombre mujer
    #>   completa     21    24
    #>   parcial      21    25

    plot(tabla1, col =c("red", "blue"), main = "jornada vs. género")

7) ¿Qué relación existe entre jornada y contrato? Análogamente, contesta a las preguntas del apartado anterior.

    tabla2 <- table(riesgos$jornada, riesgos$contrato)
    tabla2

    #>
    #>            fijo temporal
    #>   completa   38        7
    #>   parcial     0       46

    plot(tabla2, col =c("green", "blue"), main = "jornada vs. contrato")

## D. Inferencia de una proporción

8) Obtener un intervalo, con una confianza del 90 %, para el porcentaje de trabajadores con una edad menor de 30 años (2 decimales)

Definimos éxito cuando el trabajador presenta una edad < 30 años

    round(prop.test(x=25, n=91,  conf.level= 0.90)$conf.int,2)

    #> [1] 0.20 0.36
    #> attr(,"conf.level")
    #> [1] 0.9

9) En la población de trabajadores de este estudio, ¿podemos admitir que el porcentaje de individuos con una edad < 30 años es igual al 40%?. Realiza una prueba de hipótesis para responder a esta pregunta, considerando un nivel de confianza del 90%. ¿El test planteado es de una o de dos colas?
Definimos éxito cuando el trabajador presenta una edad < 30 años

El parámetro de interés es la proporción poblacional: 40%

Necesitamos saber el número de personas con menos de 35 años en nuestra muestra:

    table(edad < 30)

    #>
    #> FALSE  TRUE
    #>    66    25

    prop.test(x=25, n=91, p=0.4, alt ="two.sided", conf.level= 0.90)

    #>
    #>  1-sample proportions test with continuity correction
    #>
    #> data:  25 out of 91, null probability 0.4
    #> X-squared = 5.44, df = 1, p-value = 0.01968
    #> alternative hypothesis: true p is not equal to 0.4
    #> 90 percent confidence interval:
    #>  0.200230 0.363141
    #> sample estimates:
    #>         p
    #> 0.2747253

10) Seguimos con el apartado anterior. Algunas preguntas más:

- ¿Podemos admitir que el porcentaje de individuos con una edad < 30 años es menor del 40%?
- ¿Podemos admitir que el porcentaje de individuos con una edad < 30 años es mayor del 40%?
- ¿Los tests planteados son de una o de dos colas?

    prop.test(x=25, n=91, p=0.4, alt ="less", conf.level= 0.90)

    #>
    #>  1-sample proportions test with continuity correction
    #>
    #> data:  25 out of 91, null probability 0.4
    #> X-squared = 5.44, df = 1, p-value = 0.00984
    #> alternative hypothesis: true p is less than 0.4
    #> 90 percent confidence interval:
    #>  0.00000 0.34404
    #> sample estimates:
    #>         p
    #> 0.2747253

    prop.test(x=25, n=91, p=0.4, alt ="greater", conf.level= 0.90)

    #>
    #>  1-sample proportions test with continuity correction
    #>
    #> data:  25 out of 91, null probability 0.4
    #> X-squared = 5.44, df = 1, p-value = 0.9902
    #> alternative hypothesis: true p is greater than 0.4
    #> 90 percent confidence interval:
    #>  0.2141217 1.0000000
    #> sample estimates:
    #>         p
    #> 0.2747253

11) En la muestra, ¿cuántos trabajadores son hombres con más de 25 años?

    hjovenes <- riesgos$edad > 25  &  riesgos$genero == "hombre"
    table(hjovenes)

    #> hjovenes
    #> FALSE  TRUE
    #>    51    40

12) ¿Y cuántos trabajdores son mujeres con más de 25 años?

    mjovenes <- riesgos$edad > 25   &  riesgos$genero == "mujer"
    table(mjovenes)

    #> mjovenes
    #> FALSE  TRUE
    #>    43    48

13) A partir de una edad determinada, puede tener más impacto algunas exposiciones de riesgos en los trabajadores. Por ello nos gustaría conocer si hay diferencias de edad entre hombres y mujeres, para detectar algún colectivo que debamos proteger mejor. En la población de trabajadores de centros hospitalarios en Valencia, ¿podemos admitir que el % de hombres menores de 40 años es el mismo que el % de mujeres menores de 40 años? Utiliza una prueba de hipótesis para responder a esta pregunta. Definimos éxito cuando el trabajador tiene < 40 años

El parámetro de interés es la diferencia de dos proporciones poblacionales

Necesitamos saber el número de personas < 40 años en nuestra muestra para hombres y mujeres. Ya vimos que eran 42 y 49 respectivamente. También el número total de hombres y mujeres. Ahora ya estamos listos para utilizar la correspondiente prueba de hipótesis:

    table(riesgos$genero)

    #>
    #> hombre  mujer
    #>     42     49

    prop.test(x =c(40,48), n=c(42,49), alt ="two.sided", conf.level= 0.95)

    #> Warning in prop.test(x = c(40, 48), n = c(42, 49), alt = "two.sided",
    #> conf.level = 0.95): Chi-squared approximation may be incorrect

    #>
    #>  2-sample test for equality of proportions with continuity
    #>  correction
    #>
    #> data:  c(40, 48) out of c(42, 49)
    #> X-squared = 0.018466, df = 1, p-value = 0.8919
    #> alternative hypothesis: two.sided
    #> 95 percent confidence interval:
    #>  -0.12491923  0.07049746
    #> sample estimates:
    #>    prop 1    prop 2
    #> 0.9523810 0.9795918

14) Obtener un intervalo, con una confianza del 90 %, para la diferencia de proporciones (hombres-mujeres) de trabajadores con una edad < 40 años (2 decimales).

    round(prop.test(x =c(40,48), n=c(42,49), alt ="two.sided", conf.level= 0.95)$conf.int,2)

    #> Warning in prop.test(x = c(40, 48), n = c(42, 49), alt = "two.sided",
    #> conf.level = 0.95): Chi-squared approximation may be incorrect

    #> [1] -0.12  0.07
    #> attr(,"conf.level")
    #> [1] 0.95

## E. Inferencia de medias

15) ¿Qué relación existe entre el sueldo y el género de los trabajadores?

Primero lo veremos a nivel descriptivo: gráficamente y a través de los descriptivos de la muestra

    boxplot(riesgos$sueldo ~ riesgos$genero)

Las distribuciones de los sueldos son similares, aunque un poco menores las mujeres

También lo podemos ver con los descriptivos habituales

    summary(riesgos$sueldo[riesgos$genero=="hombre"])

    #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    #>   18000   20250   23000   22620   24000   26000

    summary(riesgos$sueldo[riesgos$genero=="mujer"])

    #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    #>   18000   20000   21000   21780   23000   26000

Ya sabemos que nos dice la muestra, pero ¿serán significativas estas diferencias? 16) Obtener un intervalo, con una confianza del 90 %, para la diferencia (mujeres-hombres) del sueldo medio de las trabajadoras y trabajadores (valores enteros)

    sueldo.hombres <- riesgos$sueldo[riesgos$genero=="hombre"]
    sueldo.mujeres <- riesgos$sueldo[riesgos$genero=="mujer"]
    round(t.test(sueldo.mujeres, sueldo.hombres, conf.level = 0.9)$conf.int,0)

    #> [1] -1630   -57
    #> attr(,"conf.level")
    #> [1] 0.9

17) Plantear una prueba no paramétrica con un nivel de confianza al 90% para ver si podemos asumir que la distribución del sueldo es el mismo en las poblaciones de trabajadores hombres y mujeres:

    wilcox.test(sueldo.mujeres, sueldo.hombres, conf.level = 0.9)

    #> Warning in wilcox.test.default(sueldo.mujeres, sueldo.hombres, conf.level =
    #> 0.9): cannot compute exact p-value with ties

    #>
    #>  Wilcoxon rank sum test with continuity correction
    #>
    #> data:  sueldo.mujeres and sueldo.hombres
    #> W = 820.5, p-value = 0.09473
    #> alternative hypothesis: true location shift is not equal to 0

18) En las poblaciones de trabajadores, ¿podemos admitir que el estrés medio es el mismo entre hombres y mujeres?

Nota: existen muchas dudas respecto a si el estrés en estos grupos sigue una distribución normal. Aplicar una prueba de hipótesis para responder a la pregunta

    estres.hombres <- riesgos$estres[riesgos$genero=="hombre"]
    estres.mujeres <- riesgos$estres[riesgos$genero=="mujer"]

En promedio los hombres tienen más estrés

    summary(estres.hombres)

    #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    #>   1.000   3.250   7.000   6.619   9.000  10.000

    summary(estres.mujeres)

    #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's
    #>   2.000   2.000   2.500   3.125   3.000   9.000       1

no hay homocedasticidad

    sd(estres.hombres)

    #> [1] 2.912985

    sd(estres.mujeres)

    #> [1] NA

Parece que algún valor perdido produce un problema en el cálculo de la SD. Lo solucionaremos así:

    sd(estres.mujeres, na.rm = TRUE)

    #> [1] 1.93099

No parece existir normalidad

    qqnorm(estres.hombres);qqline(estres.hombres)

    qqnorm(estres.mujeres);qqline(estres.mujeres)

    wilcox.test(estres.hombres, estres.mujeres)

    #> Warning in wilcox.test.default(estres.hombres, estres.mujeres): cannot
    #> compute exact p-value with ties

    #>
    #>  Wilcoxon rank sum test with continuity correction
    #>
    #> data:  estres.hombres and estres.mujeres
    #> W = 1617, p-value = 3.843e-07
    #> alternative hypothesis: true location shift is not equal to 0

## F. Inferencia para tablas de contingencia

19) Crear una nueva variable que agrupe los datos de la edad de los trabajadores en los siguientes tramos de edad: (22,30] (30,34] (34,39]:

grupo.edad = cut(riesgos$edad, breaks=c(22,30,34,39))

Obtener la tabla de frecuencias observadas para esta variable: table(grupo.edad)

    grupo.edad = cut(riesgos$edad, breaks=c(22,30,34,39))
    table(grupo.edad)

    #> grupo.edad
    #> (22,30] (30,34] (34,39]
    #>      36      37      18

20) Realizar una prueba de hipótesis para ver si podemos admitir que las variables edad (en tramos) y el sexo de los participantes son independientes. Para ello: crear la tabla de frecuencias cruzadas: tabla=table(genero, grupo.edad) y aplicar un test para tablas de contingencia.

- ¿Cuál es el valor del estadístico?
- ¿Y el valor de p del test (3 decimales)?
- ¿Cuál es tu conclusión?

    tabla=table(genero, grupo.edad)
    tabla

    #>         grupo.edad
    #> genero   (22,30] (30,34] (34,39]
    #>   hombre      15      16      11
    #>   mujer       21      21       7

    chisq.test(tabla)

    #>
    #>  Pearson's Chi-squared test
    #>
    #> data:  tabla
    #> X-squared = 2.0382, df = 2, p-value = 0.3609

    round(chisq.test(tabla)$p.value,3)

    #> [1] 0.361

21) Obtener la tabla de frecuencias que esperaríamos obtener si las variables fuesen independientes (2 decimales). ¿Es parecida a la tabla de frecuencias observadas?:

    chisq.test(table(genero, grupo.edad))$expected

    #>         grupo.edad
    #> genero    (22,30]  (30,34]  (34,39]
    #>   hombre 16.61538 17.07692 8.307692
    #>   mujer  19.38462 19.92308 9.692308

    round(chisq.test(table(genero, grupo.edad))$expected,2)

    #>         grupo.edad
    #> genero   (22,30] (30,34] (34,39]
    #>   hombre   16.62   17.08    8.31
    #>   mujer    19.38   19.92    9.69

    tabla

    #>         grupo.edad
    #> genero   (22,30] (30,34] (34,39]
    #>   hombre      15      16      11
    #>   mujer       21      21       7

## G. Bootstrap

22) Queremos obtener un intervalo de confianza del 95 % para la mediana del nivel de estrés de las trabajadores, mediante bootstrap. Para ello, seleccionaremos sólo los niveles de estrés de las trabajadores hombres: estres.hombres = riesgos$estres[riesgos$genero == “hombre”]

Creamos un vector, para almacenar la mediana de cada réplica (haremos 5000 réplicas): medianas=numeric(5000)

Fijad seed al valor 123, para obtener todos los mismos resultados: set.seed(123)

Generamos 5000 replicas de la muestra original y calculamos la mediana de cada muestra:

for(i in 1:5000) medianas[i]=median(sample(estres.hombres, replace=T))

    estres.hombres = riesgos$estres[riesgos$genero == "hombre"]
    medianas=numeric(5000)
    class(medianas)

    #> [1] "numeric"

    head(medianas)

    #> [1] 0 0 0 0 0 0

    length(medianas)

    #> [1] 5000

    set.seed(123)
    for(i in 1:5000) medianas[i]=median(sample(estres.hombres, replace=T))
    head(medianas)

    #> [1] 7 9 7 8 8 7

    summary(medianas)

    #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    #>   3.500   7.000   7.000   7.533   8.000   9.000

Determinamos un intervalo de confianza del 95% para las medianas:

    t.test(medianas)$conf.int

    #> [1] 7.50883 7.55657
    #> attr(,"conf.level")
    #> [1] 0.95

23) ¿Crees que la distribución de las medianas sigue una distribución normal? ¿Por qué?

    hist(medianas)

    boxplot(medianas)

    qqnorm(medianas);qqline(medianas)

    summary(medianas)

    #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    #>   3.500   7.000   7.000   7.533   8.000   9.000

    shapiro.test(medianas)

    #>
    #>  Shapiro-Wilk normality test
    #>
    #> data:  medianas
    #> W = 0.6575, p-value < 2.2e-16

NO es normal, luego este IC no sería adecuado 24) Finalmente obtenemos el intervalo, con una confianza del 95 %, para la mediana del estrés de los trabajadores hombres: quantile(medianas, c(0.025, 0.975)) (2 decimales)

round(quantile(medianas, c(0.025,0.975)),2) #intervalo basado en los percentiles

    round(quantile(medianas, c(0.025,0.975)),2)

    #>  2.5% 97.5%
    #>     7     9

25) Calcula el sesgo de este estadístico, frente al valor de la mediana de la muestra original.

    median(estres.hombres)

    #> [1] 7

    mean(medianas)

    #> [1] 7.5327

    sesgo = abs(median(estres.hombres) - mean(medianas))
    round(sesgo,3)

    #> [1] 0.533

26) Obtener mediante bootstrap, realizando 5000 réplicas, un intervalo de confianza del 95 % para la desviación típica del nivel de estrés de las trabajadores hombres del estudio. Fijar seed al valor 123, para obtener todos los mismos resultados (2 decimales):

set.seed(123)

    estres.hombres = riesgos$estres[riesgos$genero == "hombre"]

Creamos un vector, para almacenar la sd de cada réplica (haremos 5000 réplicas):

    desviaciones=numeric(5000)
    set.seed(123)

Generamos 5000 replicas de la muestra original y calculamos la mediana de cada muestra.

    for(i in 1:5000) desviaciones[i]=sd(sample(estres.hombres, replace=T))
    head(desviaciones)

    #> [1] 3.000290 2.297931 3.107445 2.681398 2.721276 2.724795

    summary(desviaciones)

    #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    #>   1.861   2.734   2.889   2.867   3.020   3.501

    round(quantile(desviaciones, c(0.025,0.975)),2)

    #>  2.5% 97.5%
    #>  2.39  3.25

    sesgo = abs(sd(estres.hombres) - mean(desviaciones))
    round(sesgo,3)

    #> [1] 0.046

27) ¿La distribución bootstrap de las desviaciones típicas es normal?¿Por qué?

    boxplot(desviaciones)

    hist(desviaciones)

    qqnorm(desviaciones); qqline(desviaciones)

    summary(desviaciones)

    #>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    #>   1.861   2.734   2.889   2.867   3.020   3.501

    shapiro.test(desviaciones)

    #>
    #>  Shapiro-Wilk normality test
    #>
    #> data:  desviaciones
    #> W = 0.98532, p-value < 2.2e-16

    date ()

    #> [1] "Sun Oct 18 17:35:15 2015"

    sessionInfo ()

    #> R version 3.2.1 (2015-06-18)
    #> Platform: x86_64-pc-linux-gnu (64-bit)
    #> Running under: Ubuntu precise (12.04.3 LTS)
    #>
    #> locale:
    #>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
    #>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
    #>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
    #>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
    #>  [9] LC_ADDRESS=C               LC_TELEPHONE=C
    #> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
    #>
    #> attached base packages:
    #> [1] stats     graphics  grDevices utils     datasets  methods   base
    #>
    #> other attached packages:
    #> [1] limma_3.22.7    boot_1.3-17     reshape_0.8.5   knitr_1.11
    #> [5] rmarkdown_0.8.1
    #>
    #> loaded via a namespace (and not attached):
    #>  [1] plyr_1.8.3      magrittr_1.5    markdown_0.7.7  tools_3.2.1
    #>  [5] htmltools_0.2.6 Rcpp_0.12.0     stringi_0.5-5   stringr_1.0.0
    #>  [9] digest_0.6.8    mime_0.4        evaluate_0.7.2