-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathCapitulo_04.Rmd
1514 lines (995 loc) · 71.6 KB
/
Capitulo_04.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# Revisión de estadística con R {#RER}
```{r, echo = F}
options(knitr.duplicate.label = "allow")
```
```{r, 219, child="_setup.Rmd"}
```
```{r, 220, eval=knitr::opts_knit$get("rmarkdown.pandoc.to") == "html", results='asis', echo=FALSE}
cat('<hr style="background-color:#03193b;height:2px">')
```
Esta sección revisa conceptos estadísticos importantes:
- Estimación de parámetros poblacionales desconocidos.
- Evaluación de la hipótesis.
- Intervalos de confianza.
Estos métodos se utilizan mucho en econometría. Se discutirán en el contexto simple de inferencia sobre una media poblacional desconocida, así como de diferentes aplicaciones en **R**. Dichas aplicaciones en **R** encuentran fundamento en los siguientes paquetes que no forman parte de la versión base de **R**:
+ **readxl** - permite importar datos de *Excel* a **R**.
+ **dplyr** - proporciona una gramática flexible para la manipulación de datos.
+ **MASS** - una colección de funciones para estadísticas aplicadas.
Asegúrese de que estén instalados antes de continuar e intentar replicar los ejemplos. La forma más segura de hacerlo es verificando si el siguiente fragmento de código se ejecuta sin errores.
```{r, 221, warning=FALSE, message=FALSE, eval=FALSE}
library(dplyr)
library(MASS)
library(readxl)
```
## Estimación de la media poblacional
```{r, 222, eval = my_output == "html", results='asis', echo=F, purl=F}
cat('
<div class = "keyconcept" id="KC3.1">
<h3 class = "right"> Concepto clave 3.1 </h3>
<h3 class = "left"> Estimadores y estimaciones </h3>
Los *estimadores* son funciones extraídas de una muestra de datos que parten de una población desconocida. Las *estimaciones* son valores numéricos calculados por estimadores basados en los datos de la muestra. Los estimadores son variables aleatorias porque son funciones de datos *aleatorios*. Las estimaciones son números no aleatorios.
</div>
')
```
```{r, 223, eval = my_output == "latex", results='asis', echo=F, purl=F}
cat('\\begin{keyconcepts}[Estimadores y estimaciones]{3.1}
Los \\textit{estimadores} son funciones extraídas de una muestra de datos que parten de una población desconocida. Las \\textit{estimaciones} son valores numéricos calculados por estimadores basados en los datos de la muestra. Los estimadores son variables aleatorias porque son funciones de datos \\textit{aleatorios}. Las estimaciones son números no aleatorios.
\\end{keyconcepts}')
```
Piense en alguna variable económica, por ejemplo, los ingresos por hora de los graduados universitarios, denotados por $Y$. Suponga que se está interesado en $\mu_Y$ la media de $Y$. Para calcular exactamente $\mu_Y$ se tendrían que entrevistar a todos los graduados que trabajan en el sistema económico. Simplemente no se puede hacer esto debido a limitaciones de tiempo y costos. Sin embargo, se puede extraer una muestra aleatoria de $n$ i.i.d. observaciones $Y_1, \dots, Y_n$ y estimar $\mu_Y$ usando uno de los estimadores más simples en el sentido del Concepto Clave 3.1 que uno pueda imaginar, es decir,
$$ \overline{Y} = \frac{1}{n} \sum_{i=1}^n Y_i, $$
la media muestral de $Y$. Por otra parte, se podría usar un estimador aún más simple para $\mu_Y$: la primera observación de la muestra, $Y_1$. ¿Es $Y_1$ un buen estimador? Por ahora, asuma que
$$ Y \sim \chi_{12}^2 $$
lo cual no es demasiado irracional ya que los ingresos por hora no son negativos y se espera que muchas ganancias por hora estén en un rango de $5€\,$ a $15€$. Además, es común que las distribuciones de ingresos estén sesgadas hacia la derecha, una propiedad de la distribución $\chi^2_{12}$.
```{r, 224, fig.align='center'}
# graficar la distribución chi_12^2
curve(dchisq(x, df=12),
from = 0,
to = 40,
ylab = "Densidad",
xlab = "Ganancias por hora en euros")
```
Ahora se extraerá una muestra de $n = 100$ observaciones y se tomará la primera observación $Y_1$ como una estimación de $\mu_Y$
```{r, 225, fig.align='center'}
# sembrar la semilla para la reproducibilidad
set.seed(1)
# muestra de la distribución chi_12^2, usar solo la primera observación
rsamp <- rchisq(n = 100, df = 12)
rsamp[1]
```
El estimado de $8.26$ no está muy lejos de $\mu_Y = 12$ pero es algo intuitivo que se podría hacer algo mejor: el estimador $Y_1$ descarta mucha información y su varianza es la varianza de la población:
$$ \text{Var}(Y_1) = \text{Var}(Y) = 2 \cdot 12 = 24 $$
Esto nos lleva a la siguiente pregunta: ¿Qué es un estimador *bueno* de un parámetro desconocido en primer lugar? Esta cuestión se aborda en los Conceptos clave 3.2 y 3.3.
```{r, 226, eval = my_output == "html", results='asis', echo=F, purl=F}
cat('
<div class = "keyconcept" id="KC3.2">
<h3 class = "right"> Concepto clave 3.2 </h3>
<h3 class = "left"> Sesgo, consistencia y eficiencia </h3>
Las características deseables de un estimador incluyen insesgabilidad, consistencia y eficiencia.
**Insesgabilidad:**
Si la media de la distribución muestral de algún estimador $\\hat\\mu_Y$ para la media poblacional $\\mu_Y$ es igual a $\\mu_Y$,
$$ E(\\hat\\mu_Y) = \\mu_Y, $$
el estimador es imparcial para $\\mu_Y$. El *sesgo* de $\\hat\\mu_Y$ entonces es $0$:
$$ E(\\hat\\mu_Y) - \\mu_Y = 0$$
**Consistencia:**
Se quiere que la incertidumbre del estimador $\\mu_Y$ disminuya a medida que aumenta el número de observaciones en la muestra. Más precisamente, se quiere que la probabilidad de que la estimación $\\hat\\mu_Y$ caiga dentro de un pequeño intervalo alrededor del valor real $\\mu_Y$ se acerque cada vez más a $1$ a medida que $n$ crece. Se escribe esto como
$$ \\hat\\mu_Y \\xrightarrow{p} \\mu_Y. $$
**Varianza y eficiencia:**
Se quiere que el estimador sea eficiente. Suponga que se tienen dos estimadores, $\\hat\\mu_Y$ y $\\overset{\\sim}{\\mu}_Y$ para un tamaño de muestra dado $n$ se sigue que
$$ E(\\hat\\mu_Y) = E(\\overset{\\sim}{\\mu}_Y) = \\mu_Y $$
pero
$$\\text{Var}(\\hat\\mu_Y) < \\text{Var}(\\overset{\\sim}{\\mu}_Y).$$
Entonces se prefiere usar $\\hat\\mu_Y$, ya que tiene una variación menor que $\\overset{\\sim}{\\mu}_Y$, lo que significa que $\\hat\\mu_Y$ es más *eficiente* en el uso de la información proporcionada por las observaciones en la muestra.
</div>
')
```
```{r, 227, eval = my_output == "latex", results='asis', echo=F, purl=F}
cat('\\begin{keyconcepts}[Sesgo\\comma consistencia y eficiencia]{3.2}
Las características deseables de un estimador incluyen insesgabilidad, consistencia y eficiencia.\\newline
\\textbf{Insesgabilidad::}
Si la media de la distribución muestral de algún estimador $\\hat\\mu_Y$ para la media poblacional $\\mu_Y$ es igual a $\\mu_Y$,
$$ E(\\hat\\mu_Y) = \\mu_Y, $$
el estimador es imparcial para $\\mu_Y$. El \\textit{sesgo} de $\\hat\\mu_Y$ entonces es $0$:
$$ E(\\hat\\mu_Y) - \\mu_Y = 0$$
\\textbf{Consistencia:}
Se quiere que la incertidumbre del estimador $\\mu_Y$ disminuya a medida que aumenta el número de observaciones en la muestra. Más precisamente, se quiere que la probabilidad de que la estimación $\\hat\\mu_Y$ caiga dentro de un pequeño intervalo alrededor del valor real $\\mu_Y$ se acerque cada vez más a $1$ a medida que $n$ crece. Se escribe esto como
$$ \\hat\\mu_Y \\xrightarrow{p} \\mu_Y. $$
\\textbf{Varianza y eficiencia:}
Se quiere que el estimador sea eficiente. Suponga que se tienen dos estimadores, $\\hat\\mu_Y$ y $\\overset{\\sim}{\\mu}_Y$ para un tamaño de muestra dado $n$ se sigue que
$$ E(\\hat\\mu_Y) = E(\\overset{\\sim}{\\mu}_Y) = \\mu_Y $$
pero
$$\\text{Var}(\\hat\\mu_Y) < \\text{Var}(\\overset{\\sim}{\\mu}_Y).$$
Entonces se prefiere usar $\\hat\\mu_Y$, ya que tiene una variación menor que $\\overset{\\sim}{\\mu}_Y$, lo que significa que $\\hat\\mu_Y$ es más \\textit{eficiente} en el uso de la información proporcionada por las observaciones en la muestra.
\\end{keyconcepts}')
```
## Propiedades de la media muestral {#PMM}
```{block2, consistency, type='rmdknit'}
Una forma más precisa de expresar la consistencia de un estimador $\hat\mu$ para un parámetro $\mu$ es
$$ P(|\hat{\mu} - \mu|<\epsilon) \xrightarrow[n \rightarrow \infty]{p} 1 \quad \text{for any}\quad\epsilon>0.$$
Esta expresión implica que la probabilidad de observar una desviación del valor verdadero de $\mu$ que es menor que algún $\epsilon > 0$ arbitrario converge a $1$ a medida que $n$ crece. La consistencia no requiere imparcialidad.
```
Para examinar las propiedades de la media muestral como estimador de la media poblacional correspondiente, considere el siguiente ejemplo en **R**.
Se genera una población **pop** que consta de observaciones $Y_i$, $i=1,\dots,10000$ que se originan en una distribución normal con media $\mu = 10$ y varianza $\sigma^2 = 1$.
Para investigar el comportamiento del estimador $\hat{\mu} = \bar{Y}$ se pueden extraer muestras aleatorias de esta población y calcular $\bar{Y}$ para cada una de ellas. Esto se hace fácilmente haciendo uso de la función **replicate()**. El argumento **expr** se evalúa **n** veces. En este caso, se extraen muestras de tamaños $n = 5$ y $n = 25$, se calculan las medias muestrales y se repite esto exactamente $N = 25000$ veces.
Para fines de comparación, se almacenan los resultados para el estimador $Y_1$, la primera observación en una muestra para una muestra de tamaño $5$, por separado.
```{r, 228, echo = T, eval = T, message = F, warning = F}
# generar una población ficticia
pop <- rnorm(10000, 10, 1)
# muestra de la población y estimación de la media
est1 <- replicate(expr = mean(sample(x = pop, size = 5)), n = 25000)
est2 <- replicate(expr = mean(sample(x = pop, size = 25)), n = 25000)
fo <- replicate(expr = sample(x = pop, size = 5)[1], n = 25000)
```
Comprobar que **est1** y **est2** son vectores de longitud $25000$:
```{r, 229}
# comprobar si el tipo de objeto es vector
is.vector(est1)
is.vector(est2)
# comprobar la longitud
length(est1)
length(est2)
```
El fragmento de código a continuación produce una gráfica de las distribuciones muestrales de los estimadores $\bar{Y}$ y $Y_1$ sobre la base de las muestras de $25000$ en cada caso. También se grafica la función de densidad de la distribución $\mathcal{N}(10,1)$.
```{r, 230, echo = T, eval = T, message = F, warning = F, fig.align='center'}
# graficar la densidad de la estimación Y_1
plot(density(fo),
col = "green",
lwd = 2,
ylim = c(0, 2),
xlab = "Estimados",
main = "Distribuciones muestrales de estimadores insesgados")
# agregar la estimación de densidad para la distribución de la media muestral con n = 5 a la gráfica
lines(density(est1),
col = "steelblue",
lwd = 2,
bty = "l")
# agregar la estimación de densidad para la distribución de la media muestral con n = 25 a la gráfica
lines(density(est2),
col = "red2",
lwd = 2)
# agregar una línea vertical en el parámetro verdadero
abline(v = 10, lty = 2)
# agregar la densidad N(10, 1) a la gráfica
curve(dnorm(x, mean = 10),
lwd = 2,
lty = 2,
add = T)
# agregar una leyenda
legend("topleft",
legend = c("N(10,1)",
expression(Y[1]),
expression(bar(Y) ~ n == 5),
expression(bar(Y) ~ n == 25)
),
lty = c(2, 1, 1, 1),
col = c("black","green", "steelblue", "red2"),
lwd = 2)
```
Primero, *todas* las distribuciones de muestreo (representadas por líneas continuas) se centran alrededor de $\mu = 10$. Esto es evidencia de la *imparcialidad* de $Y_1$, $\overline{Y}_{5}$ y $\overline{Y}_{25}$. Por supuesto, la densidad teórica $\mathcal{N}(10,1)$ también se centra en $10$.
A continuación, se observa la extensión de las distribuciones muestrales. Varias cosas son dignas de mención:
- La distribución de muestreo de $Y_1$ (curva verde) rastrea la densidad de la distribución $\mathcal{N}(10,1)$ (línea discontinua negra) bastante de cerca. De hecho, la distribución muestral de $Y_1$ es la distribución $\mathcal{N}(10,1)$. Esto es menos sorprendente si se tiene en cuenta que el estimador $Y_1$ no hace más que informar una observación que se selecciona al azar de una población con distribución $\mathcal{N}(10,1)$. Por lo tanto, $Y_1 \sim \mathcal{N}(10,1)$. Se debe tener en cuenta que este resultado no depende del tamaño de la muestra $n$: la distribución muestral de $Y_1$ *es siempre* la distribución de la población, sin importar el tamaño de la muestra. $Y_1$ es una buena estimación de $\mu_Y$, pero se puede hacer mejor.
- Ambas distribuciones de muestreo de $\overline{Y}$ muestran menos dispersión que la distribución de muestreo de $Y_1$. Lo cual implica que $\overline{Y}$ tiene una variación menor que $Y_1$. En vista de los Conceptos clave 3.2 y 3.3, se encuentra que $\overline{Y}$ es un estimador más eficiente que $Y_1$. De hecho, esto es válido para todos los $n > 1$.
- $\overline{Y}$ muestra un comportamiento que ilustra la coherencia (ver Concepto clave 3.2). Las densidades azul y roja están mucho más concentradas alrededor de $\mu=10$ que la verde. A medida que aumenta el número de observaciones de $1$ a $5$, la distribución muestral se ajusta en torno al parámetro verdadero. Al aumentar el tamaño de la muestra a $25$, este efecto se vuelve más evidente. Esto implica que la probabilidad de obtener estimaciones cercanas al valor real aumenta con $n$. Esto también se refleja en los valores estimados de la función de densidad cercanos a $10$: cuanto mayor es el tamaño de la muestra, mayor es el valor de la densidad.
Se recomienda que siga adelante y modifique el código. Pruebe diferentes valores para el tamaño de la muestra y vea cómo cambia la distribución muestral de $\overline{Y}$.
#### $\overline{Y}$ es el estimador de mínimos cuadrados de $\mu_Y$ {-}
Suponga que tiene algunas observaciones $Y_1,\dots,Y_n$ en $Y \sim \mathcal{N}(10,1)$ (que se desconoce) y le gustaría encontrar un estimador $m$ que prediga las observaciones también como sea posible. Por bueno se quiere decir elegir $m$ de manera que la desviación al cuadrado total entre el valor predicho y los valores observados sea pequeña. Matemáticamente, esto significa que se quiere encontrar un $m$ que minimice
\begin{equation}
\sum_{i=1}^n (Y_i - m)^2. (\#eq:sqm)
\end{equation}
Piense en $Y_i - m$ como el error cometido al predecir $Y_i$ por $m$. También se podría minimizar la suma de las desviaciones absolutas de $m$, pero minimizar la suma de las desviaciones al cuadrado es matemáticamente más conveniente (y conduce a un resultado diferente). Es por eso que el estimador que se está buscando se llama *estimador de mínimos cuadrados*. $m = \overline{Y}$, la media de la muestra, es este estimador.
Se puede mostrar esto generando una muestra aleatoria y graficando \@ref(eq:sqm) como una función de $m$.
```{r, 231, fig.align='center'}
# definir la función y vectorizarla
sqm <- function(m) {
sum((y-m)^2)
}
sqm <- Vectorize(sqm)
# extraer una muestra aleatoria y calcular la media
y <- rnorm(100, 10, 1)
mean(y)
```
```{r, 232, fig.align='center'}
# graficar la función objetivo
curve(sqm(x),
from = -50,
to = 70,
xlab = "m",
ylab = "sqm(m)")
# agregar una línea vertical en la media (y)
abline(v = mean(y),
lty = 2,
col = "darkred")
# agregar anotación en la media (y)
text(x = mean(y),
y = 0,
labels = paste(round(mean(y), 2)))
```
Observe que \@ref(eq:sqm) es una función cuadrática, por lo que solo hay un mínimo. La gráfica muestra que este mínimo se encuentra exactamente en la media muestral de los datos muestrales.
```{block2, vecfunction, type='rmdknit'}
Algunas funciones en <tt>R</tt> solo pueden interactuar con funciones que toman un vector como entrada y evalúan el cuerpo de la función en cada entrada del vector, por ejemplo <tt>curve()</tt>. A estas funciones se les da el nombre de funciones vectorizadas y, a menudo, es una buena idea escribir las funciones vectorizadas por usted mismo, aunque esto es engorroso en algunos casos. Tener una función vectorizada en <tt>R</tt> nunca es un inconveniente ya que estas funciones sirven tanto para valores únicos, así como para vectores.
Se puede ver la función <tt>sqm()</tt>, que no está vectorizada:
<tt>
sqm <- function(m) {
sum((y-m)^2) #cuerpo de la función
}
</tt>
Proporcionar, por ejemplo, <tt>c (1,2,3)</tt> como argumento <tt>m</tt> causaría un error, ya que entonces la operación <tt>ym</tt> no es válida: los vectores <tt>y</tt> y <tt>m</tt> son de dimensiones incompatibles. Es por eso que no se puede usar <tt>sqm()</tt> junto con <tt>curve()</tt>.
Aquí entra en juego <tt>Vectorize()</tt>. Genera una versión vectorizada de una función no vectorizada.
```
#### ¿Por qué es importante el muestreo aleatorio? {-}
Hasta ahora, se asume (a veces implícitamente) que los datos observados $Y_1, \dots, Y_n$ son el resultado de un proceso de muestreo que satisface el supuesto de muestreo aleatorio simple. Esta suposición a menudo se cumple al estimar la media de una población usando $\overline{Y}$. Si este no es el caso, las estimaciones pueden estar sesgadas.
Volviendo a **pop**, la población ficticia de observaciones de $10000$ y calcular la media de la población $\mu_{\texttt{pop}}$:
```{r, 233}
# calcular la media poblacional de pop
mean(pop)
```
A continuación, se toman muestras de $10$ observaciones de **pop** con **sample()** y se estima $\mu_{\texttt{pop}}$ usando $\overline{Y}$ repetidamente. Sin embargo, ahora se usa un esquema de muestreo que se desvía del muestreo aleatorio simple: en lugar de asegurar que cada miembro de la población tenga la misma probabilidad de terminar en una muestra, se asigna una probabilidad más alta de ser muestreada a las observaciones más pequeñas de $2500$ de la población estableciendo el argumento **prob** en un vector adecuado de ponderaciones de probabilidad:
```{r, 234}
# simular resultados para la media muestral cuando la suposición de variables aleatorias i.i.d. falla
est3 <- replicate(n = 25000,
expr = mean(sample(x = sort(pop),
size = 10,
prob = c(rep(4, 2500), rep(1, 7500)))))
# calcular la media muestral de los resultados
mean(est3)
```
A continuación, se traza la distribución muestral de $\overline{Y}$ para el caso de no i.i.d. y se compara con la distribución muestral cuando la suposición de i.i.d. se mantiene.
```{r, 235, fig.align='center'}
# distribución muestral de la media muestral, i.i.d. sostiene, n = 25
plot(density(est2),
col = "steelblue",
lwd = 2,
xlim = c(8, 11),
xlab = "Estimados",
main = "Cuando la suposición de i.i.d. falla")
# distribución muestral de la media muestral, i.i.d. falla, n = 25
lines(density(est3),
col = "red2",
lwd = 2)
# agrega una leyenda
legend("topleft",
legend = c(expression(bar(Y)[n == 25]~", i.i.d. falla"),
expression(bar(Y)[n == 25]~", i.i.d. se mantiene")
),
lty = c(1, 1),
col = c("red2", "steelblue"),
lwd = 2)
```
Aquí, el fracaso de la suposición de i.i.d. implica que, en promedio, se ha *subestimado* $\mu_Y$ usando $\overline{Y}$: la distribución correspondiente de $\overline{Y}$ se desplaza hacia la izquierda. En otras palabras, $\overline{Y}$ es un estimador *sesgado* para $\mu_Y$ si la suposición de i.i.d. no se sostiene.
## Pruebas de hipótesis relativas a la media de la población
En esta sección, se revisan brevemente los conceptos de la prueba de hipótesis y se discute cómo realizar pruebas de hipótesis en **R**. El objetivo en hacer inferencias sobre una media poblacional desconocida.
#### Acerca de las hipótesis y las pruebas de hipótesis {-}
En una prueba de significancia se quiere aprovechar la información contenida en una muestra como evidencia a favor o en contra de una hipótesis. Esencialmente, las hipótesis son preguntas simples que pueden responderse con un "sí" o un "no". En una prueba de hipótesis, normalmente se trata con dos hipótesis diferentes:
- La *hipótesis nula*, denotada $H_0$, es la hipótesis que se está interesado en probar.
- Debe haber una *hipótesis alternativa*, denotada $H_1$, la hipótesis que se cree que se cumple si se rechaza la hipótesis nula.
La hipótesis nula de que la media poblacional de $Y$ es igual al valor $\mu_{Y,0}$ se escribe como
$$ H_0: E(Y) = \mu_{Y,0}. $$
A menudo, la hipótesis alternativa elegida es la más general,
$$ H_1: E(Y) \neq \mu_{Y,0}, $$
lo que implica que $E(Y)$ puede ser cualquier cosa menos el valor de la hipótesis nula. Esto se llama una alternativa *de dos caras*.
En aras de la brevedad, solo se consideran alternativas de dos caras en las secciones siguientes de este apartado.
### El valor p (p-Value) {-}
Suponga que la hipótesis nula es *verdadera*. El valor $p$ es la probabilidad de extraer datos y observar un estadístico de prueba correspondiente que sea al menos tan adverso a lo que se establece bajo la hipótesis nula, como el estadístico de prueba realmente calculado usando los datos de la muestra.
En el contexto de la media poblacional y la media muestral, esta definición puede expresarse matemáticamente de la siguiente manera:
\begin{equation}
p \text{-value} = P_{H_0}\left[ \lvert \overline{Y} - \mu_{Y,0} \rvert > \lvert \overline{Y}^{act} - \mu_{Y,0} \rvert \right] (\#eq:pvalue)
\end{equation}
En \@ref(eq:pvalue), $\overline{Y}^{act}$ es la media muestral de los datos disponibles (un valor). Para calcular el valor $p$ como en \@ref(eq:pvalue), se requiere el conocimiento sobre la distribución muestral de $\overline{Y}$ (una variable aleatoria) cuando la hipótesis nula es verdadera (la *distribución nula*). Sin embargo, en la mayoría de los casos se desconoce la distribución muestral y, por tanto, la distribución nula de $\overline{Y}$. Afortunadamente, el TLC (ver Concepto clave 2.7) permite la aproximación de muestra grande
$$ \overline{Y} \approx \mathcal{N}(\mu_{Y,0}, \, \sigma^2_{\overline{Y}}) \ \ , \ \ \sigma^2_{\overline{Y}} = \frac{\sigma_Y^2}{n}, $$
asumiendo que la hipótesis nula $H_0: E(Y) = \mu_{Y, 0}$ es cierta. Con algo de álgebra se sigue para una $n$ grande que
$$ \frac{\overline{Y} - \mu_{Y,0}}{\sigma_Y/\sqrt{n}} \sim \mathcal{N}(0,1). $$
Entonces, en muestras grandes, el valor $p$ se puede calcular *sin* conocimiento de la distribución muestral exacta de $\overline{Y}$ usando la aproximación normal anterior.
### Cálculo del valor p cuando se conoce la desviación estándar {-}
Por ahora, se supone que se conoce $\sigma_{\overline{Y}}$. Entonces, se puede reescribir \@ref(eq:pvalue) como:
\begin{align}
p \text{-value} =& \, P_{H_0}\left[ \left\lvert \frac{\overline{Y} - \mu_{Y,0}}{\sigma_{\overline{Y}}} \right\rvert > \left\lvert \frac{\overline{Y}^{act} - \mu_{Y,0}}{\sigma_{\overline{Y}}} \right\rvert \right] \\
=& \, 2 \cdot \Phi \left[ - \left\lvert \frac{\overline{Y}^{act} - \mu_{Y,0}}{\sigma_{\overline{Y}}} \right\rvert\right]. (\#eq:pvaluenorm1)
\end{align}
El valor $p$ es el área en las colas de la distribución $\mathcal{N}(0,1)$ que se encuentra más allá de
\begin{equation}
\pm \left\lvert \frac{\overline{Y}^{act} - \mu_{Y,0}}{\sigma_{\overline{Y}}} \right\rvert (\#eq:pvaluenorm2)
\end{equation}
Ahora usando **R** para visualizar lo que se indica en \@ref(eq:pvaluenorm1) y \@ref(eq:pvaluenorm2). El siguiente fragmento de código replica la Figura 3.1.
```{r, 236, fig.align='center'}
# graficar la densidad normal estándar en el intervalo [-4,4]
curve(dnorm(x),
xlim = c(-4, 4),
main = "Calcular un valor p",
yaxs = "i",
xlab = "z",
ylab = "",
lwd = 2,
axes = "F")
# agregar eje x
axis(1,
at = c(-1.5, 0, 1.5),
padj = 0.75,
labels = c(expression(-frac(bar(Y)^"act"~-~bar(mu)[Y,0], sigma[bar(Y)])),
0,
expression(frac(bar(Y)^"act"~-~bar(mu)[Y,0], sigma[bar(Y)]))))
# sombrear la región del valor p/2 en la cola izquierda
polygon(x = c(-6, seq(-6, -1.5, 0.01), -1.5),
y = c(0, dnorm(seq(-6, -1.5, 0.01)),0),
col = "steelblue")
# sombrear la región del valor p/2 en la cola derecha
polygon(x = c(1.5, seq(1.5, 6, 0.01), 6),
y = c(0, dnorm(seq(1.5, 6, 0.01)), 0),
col = "steelblue")
```
### Varianza de muestra, desviación estándar de muestra y error estándar {-}
Si se desconoce $\sigma^2_Y$, debe estimarse. Esto se puede hacer usando la varianza de la muestra:
\begin{equation}
s_Y^2 = \frac{1}{n-1} \sum_{i=1}^n (Y_i - \overline{Y})^2.
\end{equation}
Además
\begin{equation}
s_Y = \sqrt{\frac{1}{n-1} \sum_{i=1}^n (Y_i - \overline{Y})^2}
\end{equation}
es un estimador adecuado para la desviación estándar de $Y$. En **R**, $s_Y$ se implementa en la función **sd()**, ver `?Sd`.
Usando **R** se puede ilustrar que $s_Y$ es un estimador consistente para $\sigma_Y$, es decir
$$ s_Y \overset{p}{\longrightarrow} \sigma_Y. $$
La idea aquí es generar una gran cantidad de muestras $Y_1,\dots,Y_n$ donde, $Y\sim \mathcal{N}(10, 9)$ digamos, estimar $\sigma_Y$ usando $s_Y$ e investigar cómo la distribución de $s_Y$ cambia a medida que $n$ aumenta.
```{r, 237, fig.align='center'}
# vector de tamaños de muestra
n <- c(10000, 5000, 2000, 1000, 500)
# observaciones de muestra, estimando usando 'sd()' y graficar las distribuciones estimadas
sq_y <- replicate(n = 10000, expr = sd(rnorm(n[1], 10, 3)))
plot(density(sq_y),
main = expression("Distribuciones de muestreo o" ~ s[Y]),
xlab = expression(s[y]),
lwd = 2)
for (i in 2:length(n)) {
sq_y <- replicate(n = 10000, expr = sd(rnorm(n[i], 10, 3)))
lines(density(sq_y),
col = i,
lwd = 2)
}
# agrega una leyenda
legend("topleft",
legend = c(expression(n == 10000),
expression(n == 5000),
expression(n == 2000),
expression(n == 1000),
expression(n == 500)),
col = 1:5,
lwd = 2)
```
El gráfico muestra que la distribución de $s_Y$ se ajusta alrededor del valor real $\sigma_Y = 3$ a medida que aumenta $n$.
La función que estima la desviación estándar de un estimador se llama *error estándar del estimador*. El concepto clave 3.4 resume la terminología en el contexto de la media muestral.
```{r, 238, eval = my_output == "html", results='asis', echo=F, purl=F}
cat('
<div class = "keyconcept" id="KC3.4">
<h3 class = "right"> Concepto clave 3.4 </h3>
<h3 class = "left"> El error estándar de $\\overline{Y}$ </h3>
Tomar una muestra i.i.d. $Y_1, \\dots, Y_n$. La media de $Y$ se estima consistentemente mediante $\\overline{Y}$, la media muestral de $ Y_i $. Dado que $\\overline{Y}$ es una variable aleatoria, tiene una distribución de muestreo con varianza $\\frac{\\sigma_Y^2}{n}$.
El error estándar de $\\overline{Y}$, denotado $SE(\\overline{Y})$ es un estimador de la desviación estándar de $\\overline{Y}$:
$$ SE(\\overline{Y}) = \\hat\\sigma_{\\overline{Y}} = \\frac{s_Y}{\\sqrt{n}} $$
El signo de intercalación (\\^) sobre $\\sigma$ indica que $\\hat\\sigma_{\\overline{Y}}$ es un estimador de $\\sigma_{\\overline{Y}}$.
</div>
')
```
```{r, 239, eval = my_output == "latex", results='asis', echo=F, purl=F}
cat('\\begin{keyconcepts}[El error estándar de $\\overline{Y}$]{3.4}
Tomar una muestra i.i.d. $Y_1, \\dots, Y_n$. La media de $Y$ se estima consistentemente mediante $\\overline{Y}$, la media muestral de $ Y_i $. Dado que $\\overline{Y}$ es una variable aleatoria, tiene una distribución de muestreo con varianza $\\frac{\\sigma_Y^2}{n}$.
El error estándar de $\\overline{Y}$, denotado $SE(\\overline{Y})$ es un estimador de la desviación estándar de $\\overline{Y}$:
$$ SE(\\overline{Y}) = \\hat\\sigma_{\\overline{Y}} = \\frac{s_Y}{\\sqrt{n}} $$
El signo de intercalación (\\^) sobre $\\sigma$ indica que $\\hat\\sigma_{\\overline{Y}}$ es un estimador de $\\sigma_{\\overline{Y}}$.
\\end{keyconcepts}
')
```
Como ejemplo para respaldar el Concepto clave 3.4, considere una muestra de $n = 100$ i.i.d. observaciones de la variable distribuida de Bernoulli $Y$ con probabilidad de éxito $p=0.1$. Entonces $E(Y)=p=0.1$ y $\text{Var}(Y)=p(1-p)$. $E(Y)$ se puede estimar mediante $\overline{Y}$, que luego tiene una variación
$$ \sigma^2_{\overline{Y}} = p(1-p)/n = 0.0009 $$
y desviación estándar
$$ \sigma_{\overline{Y}} = \sqrt{p(1-p)/n} = 0.03. $$
En este caso, el error estándar de $\overline{Y}$ puede estimarse mediante
$$ SE(\overline{Y}) = \sqrt{\overline{Y}(1-\overline{Y})/n}. $$
Se comprueba si $\overline{Y}$ y $SE(\overline{Y})$ estiman los valores verdaderos respectivos, en promedio.
```{r, 240}
# extraiga 10000 muestras de tamaño 100 y estime la media de Y y
# estimar el error estándar de la media muestral
mean_estimates <- numeric(10000)
se_estimates <- numeric(10000)
for (i in 1:10000) {
s <- sample(0:1,
size = 100,
prob = c(0.9, 0.1),
replace = T)
mean_estimates[i] <- mean(s)
se_estimates[i] <- sqrt(mean(s) * (1 - mean(s)) / 100)
}
mean(mean_estimates)
mean(se_estimates)
```
Ambos estimadores parecen no tener sesgos para los parámetros verdaderos. De hecho, esto es cierto para la media de la muestra, pero no para $SE(\overline{Y})$. Sin embargo, ambos estimadores son *consistentes* para los parámetros verdaderos.
### Cálculo del valor p cuando la desviación estándar es desconocida {-}
Cuando se desconoce $\sigma_Y$, el valor de $p$ para una prueba de hipótesis sobre $\mu_Y$ usando $\overline{Y}$ se puede calcular reemplazando $\sigma_{\overline{Y}}$ en \@ref(eq:pvaluenorm1) por el error estándar $SE(\overline{Y}) = \hat\sigma_{\overline{Y}}$. Luego,
$$ p\text{-value} = 2\cdot\Phi\left(-\left\lvert \frac{\overline{Y}^{act}-\mu_{Y,0}}{SE(\overline{Y})} \right\rvert \right). $$
Esto se hace fácilmente en **R**:
```{r, 241}
# muestra y estimación, calcula el error estándar
samplemean_act <- mean(
sample(0:1,
prob = c(0.9, 0.1),
replace = T,
size = 100))
SE_samplemean <- sqrt(samplemean_act * (1 - samplemean_act) / 100)
# hipótesis nula
mean_h0 <- 0.1
# calcular el valor p
pvalue <- 2 * pnorm(- abs(samplemean_act - mean_h0) / SE_samplemean)
pvalue
```
Más adelante en el curso, se encontrarán enfoques más convenientes para obtener estadísticos $t$ y valores $p$ usando **R**.
### La estadística t {-}
En la prueba de hipótesis, el promedio de la muestra estandarizada
\begin{equation}
t = \frac{\overline{Y} - \mu_{Y,0}}{SE(\overline{Y})} (\#eq:tstat)
\end{equation}
se llama un estadístico $t$. Dicho estadístico de $t$ juega un papel importante en la prueba de hipótesis sobre $\mu_Y$. Es un ejemplo destacado de estadística de prueba.
Implícitamente, ya se ha calculado un estadístico $t$ para $\overline{Y}$ en el fragmento de código anterior.
```{r, 242}
# calcular un estadístico t para la media de la muestra
tstatistic <- (samplemean_act - mean_h0) / SE_samplemean
tstatistic
```
Usando **R** se puede ilustrar que si $\mu_{Y,0}$ es igual al valor verdadero, es decir, si la hipótesis nula es verdadera, \@ref(eq:tstat) es aproximadamente $\mathcal{N}(0,1)$ distribuido cuando $n$ es grande.
```{r, 243}
# preparar un vector vacío para estadísticas t
tstatistics <- numeric(10000)
# establecer tamaño de muestra
n <- 300
# simular 10000 estadísticos t
for (i in 1:10000) {
s <- sample(0:1,
size = n,
prob = c(0.9, 0.1),
replace = T)
tstatistics[i] <- (mean(s)-0.1)/sqrt(var(s)/n)
}
```
En la simulación anterior, se estimó la varianza de $Y_i$ usando **var(s)**. Esto es más general que `mean(s)*(1-mean(s))`, ya que esta última requiere que los datos estén distribuidos por Bernoulli y que se sepa esto.
```{r, 244, fig.align='center'}
# graficar densidad y comparar con la densidad N(0,1)
plot(density(tstatistics),
xlab = "Estadístico t",
main = "Distribución estimada del estadístico t cuando n = 300",
lwd = 2,
xlim = c(-4, 4),
col = "steelblue")
# N(0,1) densidad (discontinua)
curve(dnorm(x),
add = T,
lty = 2,
lwd = 2)
```
A juzgar por el gráfico, la aproximación normal funciona razonablemente bien para el tamaño de muestra elegido. Esta aproximación normal ya se ha utilizado en la definición del valor $p$, ver \@ref(eq:tstat).
### Prueba de hipótesis con un nivel de significancia preespecificado {-}
```{r, 245, eval = my_output == "html", results='asis', echo=F, purl=F}
cat('
<div class = "keyconcept" id="KC3.5">
<h3 class = "right"> Concepto clave 3.5 </h3>
<h3 class = "left"> La terminología de la prueba de hipótesis </h3>
En la prueba de hipótesis, son posibles dos tipos de errores:
1. La hipótesis nula *es* rechazada aunque es cierta (error de tipo I)
2. La hipótesis nula *no se* rechaza aunque es falsa (error de tipo II)
El **nivel de significancia** de la prueba es la probabilidad de cometer un error de tipo I, que se está dispuesto a aceptar de antemano. Por ejemplo, usando un nivel de significancia preespecificado de $0.05$, se rechaza la hipótesis nula si y solo si el valor $p$ es menor que $0.05$. El nivel de significancia se elige antes de realizar la prueba.
Un procedimiento equivalente es rechazar la hipótesis nula si el estadístico de prueba observado es, en términos de valor absoluto, mayor que el **valor crítico** del estadístico de prueba. El valor crítico está determinado por el nivel de significancia elegido y define dos conjuntos de valores separados que se denominan **región de aceptación** y **región de rechazo**. La región de aceptación contiene todos los valores de la estadística de prueba para los que la prueba no rechaza, mientras que la región de rechazo contiene todos los valores para los que la prueba sí rechaza.
El **valor $p$** es la probabilidad de que, en un muestreo repetido bajo las mismas condiciones, se observe un estadístico de prueba que proporcione tanta evidencia contra la hipótesis nula como el estadístico de prueba realmente observado.
La probabilidad real de que la prueba rechace la verdadera hipótesis nula se denomina **tamaño de la prueba**. En un entorno ideal, el tamaño es igual al nivel de significancia.
La probabilidad de que la prueba rechace correctamente una hipótesis nula falsa se llama **potencia**.
</div>
')
```
```{r, 246, eval = my_output == "latex", results='asis', echo=F, purl=F}
cat('\\begin{keyconcepts}[La terminología de la prueba de hipótesis]{3.5}
En la prueba de hipótesis, son posibles dos tipos de errores:\\newline
\\begin{enumerate}
\\item La hipótesis nula \\textit{es} rechazada aunque es cierta (error de tipo I)
\\item La hipótesis nula \\textit{no se} rechaza aunque es falsa (error de tipo II)
\\end{enumerate}\\vspace{0.5cm}
El \\textit{nivel de significancia} de la prueba es la probabilidad de cometer un error de tipo I, que se está dispuesto a aceptar de antemano. Por ejemplo, usando un nivel de significancia preespecificado de $0.05$, se rechaza la hipótesis nula si y solo si el valor $p$ es menor que $0.05$. El nivel de significancia se elige antes de realizar la prueba.\\newline
Un procedimiento equivalente es rechazar la hipótesis nula si el estadístico de prueba observado es, en términos de valor absoluto, mayor que el \\textit{valor crítico} del estadístico de prueba. El valor crítico está determinado por el nivel de significancia elegido y define dos conjuntos de valores separados que se denominan \\textit{región de aceptación} y \\textit{región de rechazo}. La región de aceptación contiene todos los valores de la estadística de prueba para los que la prueba no rechaza, mientras que la región de rechazo contiene todos los valores para los que la prueba sí rechaza.\\newline
El \\textit{valor $p$} es la probabilidad de que, en un muestreo repetido bajo las mismas condiciones, se observe un estadístico de prueba que proporcione tanta evidencia contra la hipótesis nula como el estadístico de prueba realmente observado.\\newline
La probabilidad real de que la prueba rechace la verdadera hipótesis nula se denomina \\textit{tamaño de la prueba}. En un entorno ideal, el tamaño es igual al nivel de significancia.\\newline
La probabilidad de que la prueba rechace correctamente una hipótesis nula falsa se llama \\textit{potencia}.
\\end{keyconcepts}
')
```
Reconsidere el **valor p** calculado más arriba:
```{r, 247}
# comprobar si el valor p < 0.05
pvalue < 0.05
```
La condición no se cumple por lo que no se rechaza correctamente la hipótesis nula.
Cuando se trabaja con un estadístico $t$ en su lugar, es equivalente a aplicar la siguiente regla:
$$ \text{Rechazar } H_0 \text{ si } \lvert t^{act} \rvert > 1.96 $$
Se rechaza la hipótesis nula al nivel de significancia de $5\%$ si la estadística de $t$ calculada se encuentra más allá del valor crítico de 1.96 en términos de valor absoluto. $1.96$ es el cuantil $0.975$ de la distribución normal estándar.
```{r, 248}
# comprobar el valor crítico
qnorm(p = 0.975)
# compruebar si la hipótesis nula se rechaza utilizando el estadístico t calculado más arriba
abs(tstatistic) > 1.96
```
Al igual que con el valor $p$, no se puede rechazar la hipótesis nula utilizando el estadístico $t$ correspondiente. El Concepto clave 3.6 resume el procedimiento de realizar una prueba de hipótesis bilateral sobre la media poblacional $E(Y)$.
```{r, 249, eval = my_output == "html", results='asis', echo=F, purl=F}
cat('
<div class = "keyconcept" id="KC3.6">
<h3 class = "right"> Concepto clave 3.6 </h3>
<h3 class = "left"> Prueba de la hipótesis $E(Y) = \\mu_{Y,0}$ contra la alternativa $E(Y) \\neq \\mu_{Y,0}$ </h3>
1. Estimar $\\mu_{Y}$ usando $\\overline{Y}$ y calculando el $SE(\\overline{Y})$, error estándar de $\\overline{Y}$.
2. Calcular el estadístico $t$.
3. Calcular el valor $p$ y rechazar la hipótesis nula en el nivel de significancia $5\\%$ si el valor $p$ es menor que $0.05$ o, de manera equivalente, si
$$ \\left\\lvert t^{act} \\right\\rvert > 1.96. $$
</div>
')
```
```{r, 250, eval = my_output == "latex", results='asis', echo=F, purl=F}
cat('\\begin{keyconcepts}[Probar la hipótesis $E(Y) = \\mu_{Y,0}$ contra la alternativa $E(Y) \\neq \\mu_{Y,0}$]{3.6}
\\begin{enumerate}
\\item Estimar $\\mu_{Y}$ usando $\\overline{Y}$ y calculando el $SE(\\overline{Y})$, error estándar de $\\overline{Y}$.
\\item Calcular el estadístico $t$.
\\item Calcular el valor $p$ y rechazar la hipótesis nula en el nivel de significancia $5\\%$ si el valor $p$ es menor que $0.05$ o, de manera equivalente, si
$$ \\left\\lvert t^{act} \\right\\rvert > 1.96. $$
\\end{enumerate}
\\end{keyconcepts}
')
```
### Alternativas unilaterales {-}
A veces interesa probar si la media es mayor o menor que algún valor hipotético bajo el nulo. Para ceñirse al curso, se toma la presunta brecha salarial entre los trabajadores con buena educación y los menos educados. Dado que se anticipó que existe tal diferencial, una alternativa relevante (a la hipótesis nula de que no existe un diferencial salarial) es que los individuos bien educados ganan más; es decir, que el salario promedio por hora para este grupo, $\mu_Y$ es *mayor* que $\mu_{Y, 0}$, el salario promedio de los trabajadores con menos educación que se asume que se conoce aquí por simplicidad (la sección \@ref(CMDP) analiza cómo probar la equivalencia de medias poblacionales desconocidas).
Este es un ejemplo de una *prueba del lado derecho* y el par de hipótesis se elige para ser
$$ H_0: \mu_Y = \mu_{Y,0} \ \ \text{vs} \ \ H_1: \mu_Y > \mu_{Y,0}. $$
Se rechaza la hipótesis nula si el estadístico de prueba calculado es mayor que el valor crítico $1.64$, el equivalente de $0.95$ de la distribución $\mathcal{N}(0,1)$. Esto asegura que $1-0.95=5\%$ masa de probabilidad permanece en el área a la derecha del valor crítico. Como antes, se puede visualizar esto en **R** usando la función **polygon()**.
```{r, 251, fig.align='center'}
# graficar la densidad normal estándar en el dominio [-4,4]
curve(dnorm(x),
xlim = c(-4, 4),
main = "Región de rechazo de una prueba del lado derecho",
yaxs = "i",
xlab = "Estadístico t",
ylab = "",
lwd = 2,
axes = "F")
# Agregar el eje x
axis(1,
at = c(-4, 0, 1.64, 4),
padj = 0.5,
labels = c("", 0, expression(Phi^-1~(.95)==1.64), ""))
# Sombrear la región de rechazo en la cola izquierda
polygon(x = c(1.64, seq(1.64, 4, 0.01), 4),
y = c(0, dnorm(seq(1.64, 4, 0.01)), 0),
col = "darkred")
```
De manera análoga, para la prueba del lado izquierdo se tiene
$$ H_0: \ mu_Y = \ mu_ {Y, 0} \ \ \ text {vs.} \ \ H_1: \ mu_Y <\ mu_ {Y, 0}. $$
El nulo se rechaza si el estadístico de prueba observado no alcanza el valor crítico que, para una prueba en el nivel de significancia de $0.05$, está dado por $-1.64$, el equivalente de $0.05$ -cuantil de la distribución $\mathcal{N}(0,1)$. $5\%$ masa de probabilidad se encuentra a la izquierda del valor crítico.
Es sencillo adaptar el fragmento de código anterior al caso de una prueba del lado izquierdo. Solo se tiene que ajustar el sombreado de color y las marcas de graduación.
```{r, 252, fig.align='center'}
# graficar la densidad normal estándar en el dominio [-4,4]
curve(dnorm(x),
xlim = c(-4, 4),
main = "Región de rechazo de una prueba del lado izquierdo",
yaxs = "i",
xlab = "Estadístico t",
ylab = "",
lwd = 2,
axes = "F")
# Agregar eje x
axis(1,
at = c(-4, 0, -1.64, 4),
padj = 0.5,
labels = c("", 0, expression(Phi^-1~(.05)==-1.64), ""))
# Región de rechazo de sombra en la cola derecha
polygon(x = c(-4, seq(-4, -1.64, 0.01), -1.64),
y = c(0, dnorm(seq(-4, -1.64, 0.01)), 0),
col = "darkred")
```
## Intervalos de confianza para la media de la población
Como se enfatizó anteriormente, nunca se calcula el valor *exacto* de la media poblacional de $Y$ usando una muestra aleatoria. Sin embargo, se pueden calcular los intervalos de confianza para la media de la población. En general, un intervalo de confianza para un parámetro desconocido es una receta que, en muestras repetidas, produce intervalos que contienen el parámetro verdadero con una probabilidad preespecificada, el *nivel de confianza*. Los intervalos de confianza se calculan utilizando la información disponible en la muestra. Dado que esta información es el resultado de un proceso aleatorio, los intervalos de confianza son variables aleatorias en sí mismas.
El Concepto clave 3.7 muestra cómo calcular los intervalos de confianza para la media poblacional desconocida $E(Y)$.
```{r, 253, eval = my_output == "html", results='asis', echo=F, purl=F}
cat('
<div class = "keyconcept" id="KC3.7">
<h3 class = "right"> Concepto clave 3.7 </h3>
<h3 class = "left"> Intervalos de confianza para la media poblacional </h3>
Un intervalo de confianza de $95\\%$ para $\\mu_Y$ es una *variable aleatoria* que contiene el verdadero $\\mu_Y$ en $95\\%$ de todas las muestras aleatorias posibles. Cuando $n$ es grande, se puede usar la aproximación normal. Entonces, $99\\%$, $95\\%$, $90\\%$ los intervalos de confianza son:
\\begin{align}
&99\\%\\text{ intervalo de confianza para } \\mu_Y = \\left[ \\overline{Y} \\pm 2.58 \\times SE(\\overline{Y}) \\right], \\\\
&95\\%\\text{ intervalo de confianza para } \\mu_Y = \\left[\\overline{Y} \\pm 1.96 \\times SE(\\overline{Y}) \\right], \\\\
&90\\%\\text{ intervalo de confianza para } \\mu_Y = \\left[ \\overline{Y} \\pm 1.64 \\times SE(\\overline{Y}) \\right].
\\end{align}
Estos intervalos de confianza son conjuntos de hipótesis nulas que no se pueden rechazar en una prueba de hipótesis bilateral con el nivel de confianza dado.
Ahora considerar las siguientes declaraciones.
1. En muestreo repetido, el intervalo
$$ \\left[ \\overline{Y} \\pm 1.96 \\times SE(\\overline{Y}) \\right] $$
cubre el valor real de $\\mu_Y$ con una probabilidad de $95\\%$.
2. Se ha calculado $\\overline{Y} = 5.1$ y $SE(\\overline{Y})=2.5$ por lo que el intervalo
$$ \\left[ 5.1 \\pm 1.96 \\times 2.5 \\right] = \\left[0.2,10\\right] $$
cubre el valor real de $\\mu_Y$ con una probabilidad de $95\\%$.
Si bien 1. es correcto (esto está en línea con la definición anterior), 2. está incorrecto y ninguno de profesionista quiere leer una oración de este tipo en un trabajo final, examen escrito o similar, creealo.
La diferencia es que, mientras que 1. es la definición de una variable aleatoria, 2. es un posible *resultado* de dicha variable aleatoria, por lo que no tiene sentido hacer ninguna afirmación probabilística al respecto. ¡O el intervalo calculado cubre $\\mu_Y$ *o no*!
</div>
')
```
```{r, 254, eval = my_output == "latex", results='asis', echo=F, purl=F}
cat('\\begin{keyconcepts}[Intervalos de confianza para la media poblacional]{3.7}
Un intervalo de confianza de $95\\%$ para $\\mu_Y$ es una \\texttt{variable aleatoria} que contiene el verdadero $\\mu_Y$ en $95\\%$ de todas las muestras aleatorias posibles. Cuando $n$ es grande, se puede usar la aproximación normal. Entonces, $99\\%$, $95\\%$, $90\\%$ los intervalos de confianza son:\\newline
\\begin{align}
&99\\%\\text{ intervalo de confianza para } \\mu_Y = \\left[ \\overline{Y} \\pm 2.58 \\times SE(\\overline{Y}) \\right], \\\\
&95\\%\\text{ intervalo de confianza para } \\mu_Y = \\left[\\overline{Y} \\pm 1.96 \\times SE(\\overline{Y}) \\right], \\\\
&90\\%\\text{ intervalo de confianza para } \\mu_Y = \\left[ \\overline{Y} \\pm 1.64 \\times SE(\\overline{Y}) \\right].
\\end{align}
Estos intervalos de confianza son conjuntos de hipótesis nulas que no se pueden rechazar en una prueba de hipótesis bilateral con el nivel de confianza dado.\\newline
Ahora considerar las siguientes declaraciones.\\newline
\\begin{enumerate}
\\item En muestreo repetido, el intervalo
$$ \\left[ \\overline{Y} \\pm 1.96 \\times SE(\\overline{Y}) \\right] $$
cubre el valor real de $\\mu_Y$ con una probabilidad de $95\\%$.
\\item Se ha calculado $\\overline{Y} = 5.1$ y $SE(\\overline{Y})=2.5$ por lo que el intervalo
$$ \\left[5.1 \\pm 1.96 \\times 2.5 \\right] = \\left[0.2,10\\right] $$
cubre el valor real de $\\mu_Y$ con una probabilidad de $95\\%$.
\\end{enumerate}\\vspace{0.5cm}
Si bien 1. es correcto (esto está en línea con la definición anterior), 2. está incorrecto y ninguno de profesionista quiere leer una oración de este tipo en un trabajo final, examen escrito o similar, creealo.
La diferencia es que, mientras que 1. es la definición de una variable aleatoria, 2. es un posible \\textit{resultado} de dicha variable aleatoria, por lo que no tiene sentido hacer ninguna afirmación probabilística al respecto. ¡O el intervalo calculado cubre $\\mu_Y$ \\textit{o no}!
\\end{keyconcepts}
')
```
En **R**, probar hipótesis sobre la media de una población sobre la base de una muestra aleatoria es muy fácil debido a funciones como **t.test()** del paquete **stats**. Produce un objeto del tipo **list**. Afortunadamente, una de las formas más sencillas de usar **t.test()** es cuando se desea obtener un intervalo de confianza de $95\%$ para alguna media poblacional. Comenzando por generar algunos datos aleatorios y llamando a **t.test()** junto con **ls()** para obtener un desglose de los componentes de salida.
```{r, 255}
# sembrar semilla
set.seed(1)
# generar algunos datos de muestra
sampledata <- rnorm(100, 10, 10)
# comprobar el tipo de resultado producido por t.test
typeof(t.test(sampledata))
# mostrar los elementos de la lista producidos por t.test
ls(t.test(sampledata))
```
Aunque se seinforman muchos elementos, por el momento solo interesa calcular un conjunto de confianza de $95\%$ para la media.
```{r, 256}
t.test(sampledata)$"conf.int"
```
Esto indica que el intervalo de confianza de $95\%$ es
$$ \left[9.31, 12.87\right]. $$
En este ejemplo, el intervalo calculado obviamente cubre el verdadero $\mu_Y$ que se sabe que es $10$.
Resulta importante echar un vistazo a toda la salida estándar producida por **t.test()**.
```{r, 257}
t.test(sampledata)
```
Se puede ver que **t.test()** no solo calcula un intervalo de confianza de $95\%$ sino que automáticamente realiza una prueba de significancia bilateral de la hipótesis $H_0: \mu_Y = 0$ al nivel de $5\%$ e informa los parámetros relevantes de la misma: la hipótesis alternativa, la media estimada, el estadístico $t$ resultante, los grados de libertad de la distribución $t$ subyacente (**t.test()** utiliza realizar la aproximación normal) y el valor $p$ correspondiente. ¡Esto es muy conveniente!
En este ejemplo, se llegó a la conclusión de que la media poblacional *es significativamente* diferente de $0$ (lo cual es correcto) al nivel de $5\%$, ya que $\mu_Y = 0$ no es un elemento del $95\%$ intervalo de confianza
$$ 0 \not\in \left[9.31,12.87\right]. $$
Se llega a un resultado equivalente cuando se usa la regla de rechazo del valor $p$, ya que
$$ p\text{-value} = 2.2\cdot 10^{-16} \ll 0.05. $$
## Comparación de medias de diferentes poblaciones {#CMDP}
Suponga que está interesado en las medias de dos poblaciones diferentes, denotadas como $\mu_1$ y $\mu_2$. Más específicamente, interesa si estas medias poblacionales son diferentes entre sí y planear usar una prueba de hipótesis para verificar esto sobre la base de datos de muestra independientes de ambas poblaciones. Un par de hipótesis adecuadas es
\begin{equation}
H_0: \mu_1 - \mu_2 = d_0 \ \ \text{vs.} \ \ H_1: \mu_1 - \mu_2 \neq d_0 (\#eq:hypmeans)
\end{equation}
donde $d_0$ denota la diferencia hipotética de medias (entonces $d_0 = 0$ cuando las medias son iguales, bajo la hipótesis nula). En el curso se enseña que $H_0$ se puede probar con el estadístico $t$
\begin{equation}
t=\frac{(\overline{Y}_1 - \overline{Y}_2) - d_0}{SE(\overline{Y}_1 - \overline{Y}_2)} (\#eq:tstatmeans)
\end{equation}
donde
\begin{equation}
SE(\overline{Y}_1 - \overline{Y}_2) = \sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}.
\end{equation}
Esto se denomina prueba $t$ de dos muestras. Para $n_1$ y $n_2$ grandes, \@ref(eq:tstatmeans) es normal estándar bajo la hipótesis nula. De manera análoga a la prueba simple $t$, se pueden calcular intervalos de confianza para la verdadera diferencia en las medias poblacionales:
$$ (\overline{Y}_1 - \overline{Y}_2) \pm 1.96 \times SE(\overline{Y}_1 - \overline{Y}_2) $$
es un intervalo de confianza de $95\%$ para $d$. <br> En **R**, las hipótesis como en \@ref(eq:hypmeans) también se pueden probar con **t.test()**. Se debe tener en cuenta que **t.test()** elige $d_0 = 0$ de forma predeterminada. En consecuenci, esto se puede cambiar configurando el argumento **mu**.
El siguiente fragmento de código demuestra cómo realizar una prueba $t$ de dos muestras en **R** utilizando datos simulados.