forked from tereom/fundamentos
-
Notifications
You must be signed in to change notification settings - Fork 0
/
03-tipos-de-estudio.Rmd
792 lines (627 loc) · 35.2 KB
/
03-tipos-de-estudio.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
# Tipos de estudio y experimentos
```{r setup, include=FALSE, message=FALSE}
library(tidyverse)
source("R/funciones_auxiliares.R")
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning=FALSE)
comma <- function(x) format(x, digits = 2, big.mark = ",")
theme_set(theme_minimal())
```
<!-- ```{block, type='ejercicio'} -->
<!-- **Pregunta de entrevista de Google [@Chihara]** -->
<!-- Imagina que eres consultor y te preguntan lo siguiente (ver siguiente figura): -->
<!-- Estoy haciendo una comparación de antes y después donde la hipótesis alternativa -->
<!-- es pre.media.error > post.media.error. La distribución de ambas muestras es -->
<!-- sesgada a la derecha. ¿Qué prueba me recomiendas para ésta situación? -->
<!-- ``` -->
<!-- ```{r grafica-pcr, echo=FALSE, warning=FALSE, message=FALSE, fig.cap = "Error CPR, gráfica de densidad.", fig.height = 3.2, fig.width = 4.2} -->
<!-- library(tidyverse) -->
<!-- theme_set(theme_minimal()) -->
<!-- pre <- tibble(group = "Pre", y = rgamma(10000, 2, 0.8)-0.2) -->
<!-- post <- tibble(group = "Post", y = rgamma(10000, 1.7, 0.8)-0.2) -->
<!-- ggplot(bind_rows(pre, post), aes(x = y, color = group)) + -->
<!-- geom_density() + -->
<!-- xlab("CPR error") + labs(color = "") -->
<!-- ``` -->
La siguiente imagen de [Roger Peng](https://simplystatistics.org/2019/04/17/tukey-design-thinking-and-better-questions/)
representa una situación común a la que se enfrenta el analista de datos, y se
desarrolló en el contexto de preguntas vagas. En el esquema hay tres caminos:
uno es ideal que pocas veces sucede,
otro produce respuestas poco útiles pero es fácil, y otro es tortuoso pero que
caracteriza el mejor trabajo de análisis de datos:
```{r, echo = FALSE, message = FALSE, fig.cap = "Adaptado de R. Peng: [Tukey, design thinking and better questions.](https://simplystatistics.org/2019/04/17/tukey-design-thinking-and-better-questions/)", warning=FALSE}
library(tidyverse)
theme_set(theme_minimal())
puntos <- tibble(x = c(0.5, 1.2, 4, 4), y = c(0.5, 4, 0.5, 5),
etiqueta = c("dónde\ncomenzamos\nrealmente", "Análisis de datos \n poco útil, de bajo impacto",
"dónde creeemos\nque comenzamos", "Nuestra \n Meta "))
set.seed(211)
browniano <- tibble(x = 0.5 + cumsum(c(0,rnorm(50, 0.03, 0.1))) ,
y = 0.5 + cumsum(c(0, rnorm(50, 0.02, 0.2))))
puntos <- bind_rows(puntos, tail(browniano, 1) %>% mutate(etiqueta = "Terminamos?!?"))
flechas <- tibble(x = c(0.5, 4), y = c(0.5, 0.5), xend = c(1.2, 4), yend = c(4, 5))
ggplot(puntos, aes(x = x, y = y)) +
xlab("Calidad de la regunta") +
ylab("Peso de la evidencia") +
theme(axis.text.x=element_blank(),
axis.text.y=element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
geom_segment(data = flechas, aes(xend=xend, yend=yend),
arrow = arrow(length = unit(0.3, "inches"))) +
geom_path(data = browniano) +
geom_point(data = browniano) +
geom_point(colour="red", size = 5) +
geom_text(aes(label = etiqueta), vjust = -0.5, hjust = 1.1, size = 4.2) +
#labs(caption = "Adaptado de R. Peng: Tukey, design thinking and better questions.") +
xlim(c(-0.1 , 4)) + ylim(c(0,6))
```
## De datos y poblaciones {-}
Los datos no son el fin último de un estudio. Son el mecanismo que podemos utilizar para poder contestar
preguntas acerca de la población que no vemos. Pensemos en la encuesta realizada en el Reino Unido
sobre parejas sexuales del sexo opuesto
que una persona en el rango de edad de 35-44 años declara tener. Este estudio está tomado de @spiegelhalter2019art.
Los datos están reportados en la encuesta Natsal-3 que puede encontrarse en [C.H. Mercer et al., ‘Changes in Sexual Attitudes and Lifestyles in Britain through the Life Course and Over Time: Findings from the National Surveys of Sexual Attitudes and Lifestyles (Natsal)’, 2013](https://www.thelancet.com/journals/lancet/article/PIIS0140-6736(13)62035-8/fulltext).
Estos datos corresponden a un total de 796 hombres y 1,193 mujeres encuestadas, y están ponderados por el diseño
estratificado de la encuesta.
```{r, echo = FALSE}
library(ggplot2)
nastal <-read.csv("data/nastal.csv", header=TRUE)
nastal <- nastal %>%
rename(NumParejas = NumPartners, ConteoH = MenCount, ConteoM = WomenCount)
sample_n(nastal, 10)
```
Como en los ejemplos anteriores, podemos calcular un resumen rápido para hombres:
```{r, echo = FALSE}
nastal %>%
summarise(Datos = rep(NumParejas,ConteoH)) %>%
pull(Datos) %>%
summary()
```
así como un resumen rápido para las mujeres encuestadas:
```{r, echo = FALSE}
nastal %>%
summarise(Datos = rep(NumParejas,ConteoM)) %>%
pull(Datos) %>%
summary()
```
Un gráfico sencillo nos ayudará a ilustrar la distirbución de las respuestas por género:
```{r, fig.align ='center', out.width='95%', fig.height = 3, fig.width = 5}
nastal %>%
select(NumParejas, ConteoH, ConteoM) %>%
rename(Hombres = ConteoH, Mujeres = ConteoM) %>%
gather('Hombres', 'Mujeres', key='Género', value = 'Conteo') %>%
ggplot(aes(x = NumParejas)) +
geom_bar(aes(y = Conteo, fill = Género), stat = 'identity', position = 'dodge') +
scale_x_continuous(breaks = c(0,5,10,15,20,25,30,35,40,45,50), limits=c(0,50)) +
scale_colour_brewer(palette = "Set1") +
labs(x="Número reportado de parejas sexuales del genero opuesto")
```
¿Cómo podemos generalizar a la población después de haber observado dichos resultados en la encuesta?
```{block, type = 'comentario'}
Podemos seguir el siguiente tren de pensamiento:
* El *número registrado* de parejas sexuales de los participantes nos habla acerca de ...
* El *número real* de parejas sexuales en nuestra *muestra*, lo que nos habla acerca de ...
* El número de parejas en la *población de estudio*, lo que nos habla acerca de ...
* El número de parejas sexuales en el Reino Unido, lo cual es la *población objetivo*.
```
```{block, type = 'comentario'}
Los puntos más débiles en la generalización son los siguientes:
* ¿Podemos asumir que los encuestados responderán de manera exacta la pregunta? Observa los *picos* en el eje horizontal.
* ¿Podemos esperar que los encuestados hayan sido escogidos de manera aleatoria de aquellos que son elegibles? Posiblemente, pero ¿podemos esperar los que aceptaron la encuesta son representativos?
* ¿Podemos asegurar que la muestra de encuestados representa la población adulta del país?
```
## Hacia el trabajo como **Juez** {-}
Este papel lo tomamos cuando queremos describir algo más allá de los datos que observamos.
Para esto necesitamos realizar *inferencia inductiva.* El peligro es que la inducción es un proceso
generalmente lleno con incertidumbre, pues implica tomar instancias muy particulares para poder emitir
juicios generales. En cambio, el trabajo *deductivo* considera una secuencia de
implicaciones lógicas que nos llevan de generalidades a casos particulares.
Posiblemente, a lo largo de su preparación, o en cursos anteriores de estadística, han
considerado los casos cuando los datos que observamos son seleccionados al azar de la población
objetivo. Sin embargo, esto raramente sucede en la vida real, y es por esto que es necesario considerar
el proceso desde la captura de datos hasta la población objetivo. Por ejemplo, la población de adultos
en Reino Unido con una vida sexual activa.
En general nos interesa que nuestros datos sean:
- Confiables, es decir, con poca variabilidad y que sea un evento repetible.
- Válidos, en el sentido que en verdad se esté midiendo lo que queremos medir, y que no haya sesgos.
Por otro lado, para poder asegurar que la muestra sea adecuada y nos permita observar
de manera fiable la población necesitamos que el estudio tenga *validez interna.* La forma
más efectiva de reducir el sesgo es por medio de **muestreo aleatorio.**
Por último, nos interesa que haya *validez externa* en los datos, lo cual significa que en verdad
nuestras unidades de observación representen la poblacion de interés.
#### Proceso generador de datos {-}
Es por esto que entre las preguntas que se debe hacer el analista de datos una fundamental es
en entender el **proceso generador de datos**, pues esto determinará qué
otras preguntas son relevantes, tanto en términos prácticos como estadísticos.
* La **inferencia estadística** busca hacer afirmaciones, cuantificadas de
manera probabilista, acerca de datos que no tenemos, usando regularidades y
conocimiento de datos que sí tenemos disponibles y métodos cuantitativos.
* Para hacer afirmaciones inferenciales **eficientes y bien calibradas** (con
garantías estadísticas de calibración) a preguntas donde queremos generalizar de
muestra a población, se requiere conocer con precisión el proceso que genera los
datos muestrales.
* Esto incluye saber con detalle cómo se seleccionaron los datos a partir de
los que se quiere hacer inferencia.
En este caso, eficiente quiere decir que aprovechamos toda la información que
está en los datos observados de manera que nuestros rangos de incertidumbre son
lo más chico posibles (además de estar correctamente calibrados).
Por su parte, probabilísticamente bien calibrados se refiere a que, lo que
decimos que puede ocurrir con 10% de probabilidad ocurre efectivamente 1 de cada
10 veces, si decimos 20% entonces ocurre 2 de 10, etc.
Veremos que para muestras dadas naturalmente, a veces es muy difiícil entender a
fondo el proceso generación de la muestra.
#### Ejemplo: Prevalencia de anemia {-}
```{r simular-datos-anemia, echo = FALSE}
# 5-9: 11,063,000
# 10-14: 11,143,000
paciente <- tibble(edad = sample(6:15, 5000, prob = c(2, 2, 2, 1, 1, 1, 0.5, 0.5,
0.2, 0.2), replace = TRUE),
padecimiento = sample(c("infección respiratoria",
"infección intestinal", "asma", "úlcera",
"picadura alacrán", "mordedura de perro",
"apendcitis"), 5000, replace = TRUE),
sexo = sample(c("hombre", "mujer"), 5000, replace = TRUE)
)
paciente <- paciente %>%
mutate(z = 1 - 0.15 * edad + case_when(padecimiento == "infección respiratoria" ~ 0.3,
padecimiento == "infección intestinal" ~ 0.7,
padecimiento == "infección respiratoria" ~ 0.2,
padecimiento == "infección respiratoria" ~ 0.9,
TRUE ~ 0),
anemia = rbinom(5000, 1, prob = exp(z) / (1 + exp(z)))
) %>%
select(-z)
```
Supongamos que nos interesa conocer el porcentaje de menores en edad escolar,
(entre 6 y 15 años), con
anemia en México. La fuente de datos disponible corresponde a registros del IMSS
de hospitalizaciones de menores, ya sea por anemia o
que por otra causa (infecciones gastrointestinales, apendicitis, tratamiento de
leucemia, ...), se registró
si el menor tenía anemia. En nuestra muestra el 47% de los niños tiene anemia.
```{r}
head(paciente)
```
- ¿Qué nos dice esta cantidad acerca de la anemia en la población de menores de edad en la república mexicana?
- ¿Podemos hacer inferencia estadística?
- ¿Cómo calculamos intervalos de confianza?
```{r}
# Si calculo el error estándar de la p estimada como sigue, es correcto?
p <- mean(paciente$anemia)
sqrt(p * (1 - p) / 5000)
```
En la situación ideal diseñaríamos una muestra aleatoria de menores de edad,
por ejemplo, utilizando el registro en educación primaria de la SEP, y
mediríamos la prevalencia de anemia en la muestra, usaríamos esta muestra para
estimar la prevalencia en la población y tendríamos además las herramientas
para medir la incertidumbre de nuestra estimación (reportar intervalos,
o errores estándar).
En el caso de prevalencia de anemia, discutiendo con médicos e investigadores
nos informan que la anemia se presenta en tasas más altas en niños más chicos.
```{r}
paciente %>%
count(edad) %>%
mutate(prop = round(100 * n / sum(n)))
```
Y consultando con las proyecciones de población notamos que los niños chicos
están sobre-representados en la muestra. Lo que nos hace considerar que debemos
buscar una manera de ponderar nuestras observaciones para que reflejen a la
población de niños en el país.
Más aún, investigamos que algunas enfermedades están asociadas a mayor
prevalencia de anemia:
```{r}
paciente %>%
count(padecimiento) %>%
arrange(-n)
```
Utilizamos esta información para modelar y *corregir* nuestra estimación
original. Por ejemplo con modelos de regresión. Sin embargo,
debemos preguntarnos:
- ¿Hay más variables qué nos falta considerar?
- Nuestras estimaciones están bien calibradas?
#### Población {-}
Hasta ahora hemos hablado de muestras de datos. Es un caso muy común en las encuestas. Sin embargo,
también es posible econtrar casos dónde tengamos acceso a todo el conjunto de datos de interés. Ejemplos
de esto son los casos de donde los casos se registran de manera continua como estudios de compras en línea, o
históricos transaccionales en un banco.
Aún en estas situaciones hay que considerar evaluar si en verdad todo lo que nos interesa se registra.
Por ejemplo, las carpetas de investigación de crímenes en la ciudad de México: ¿contienen el reporte de
todos los posibles crímenes en la ciudad?
#### Distribuciones {-}
Hasta ahora hemos mencionado el concepto de distribución como el patrón que presentan los datos
(valores centrales, dispersión, rango, etc.). A esta distribución le llamamos **distribución muestral**
o **empírica**. En general, esperamos que nuestros datos tengan las mismas características (estadísticas)
que la población de donde provienen. Por ejemplo, cuando un fenómeno es generado por pequeñas influencias
hablamos de la distribución Normal o Gaussiana a nivel teórico. El siguiente ejemplo es tomado del libro de
@spiegelhalter2019art, y es sobre el peso de los bebés al nacer para poblaciones caucásicas.
```{r, echo = FALSE, fig.align = 'center', out.width = '95%', message = FALSE}
weights=c(1500, 2000, 2500, 3000, 3500, 4000, 4500,5000)
mids=weights+250
n=c(5+48+308+1130, 12679, 124209, 442891, 389275, 108886, 14936,1345) # numbers in each bin
N=sum(n) # total number of babies
area=N*500 # number * binwidth = total area of histogram
lbw = sum(n[1:2]) # number with low birth weight (less than 2500)
lbw.percent=100*lbw/N # % low birth weight
# 1.3%
#calculate mean and sd of population
# could use sheppard's correction
birth.mean=sum(n*mids/N)
birth.sd=sqrt( sum(n*(mids-birth.mean)^2)/N)
# per cent less than 2500 from normal approximation
lbw.est = 100 * pnorm(2500,birth.mean, birth.sd)
# 1.7%, good approxmation
#25th and 75th percentiles of population
# qnorm(0.25, birth.mean,birth.sd)
# qnorm(0.75, birth.mean,birth.sd)
# percentile of baby weighing 2910
xw = 2910
# pnorm(xw, birth.mean,birth.sd)
par(mfrow=c(2,2))
# setup plot ranges noting max of normal density is at mean
xrange <- c(1500,5500)
yrange <- range( c(n, area*dnorm(birth.mean, birth.mean, birth.sd), 0))
scale=0.6
par(mar=c(5,0,1,0)+0.1)
# (a) empirical distribution and fitted normal
plot(xrange, yrange, type = "n", xlab = "", ylab = "",
bty="n",axes=F,main="(a) Distribution of birthweights", cex=scale)
axis(1,cex=scale)
# draw bars using rect and density using curve
rect(weights, 0, weights + 500, n, col = "lightblue")
curve(area*dnorm(x, birth.mean, birth.sd), min(xrange), max(xrange), add = TRUE,
lwd=3, col="blue")
lines(c(xw,xw),yrange,col="red",lwd=2)
# (b) plot with sds
plot(xrange, yrange, type = "n", xlab = "", ylab = "",
bty="n",axes=F,,main="(b) Mean +/- 1, 2, 3 SDs" )
axis(1)
curve(area*dnorm(x, birth.mean, birth.sd), min(xrange), max(xrange), add = TRUE, lwd=3, col="blue")
I=-3:3
x1=birth.mean+I*birth.sd
y1=area*dnorm(x1,birth.mean, birth.sd)
label=c("-3 SDs", "-2 SDs", "-1 SD", "mean", "+1 SD","+2 SDs", "+3 SDs")
bit=10000
xx=250
shift=c(-xx,-xx,-xx,0,xx,xx,xx)
for(i in 1:7){
lines(c(x1[i],x1[i]), c(0,y1[i]),lwd=2)
text(x1[i]+shift[i],y1[i]+bit,label[i],cex=0.75)
}
lines(c(xw,xw),yrange,col="red",lwd=2)
# (c) Percentiles
plot(xrange, yrange, type = "n", xlab = "Birthweight (gms)", ylab = "",
bty="n",axes=F,,main="(c) Percentiles" )
axis(1)
curve(area*dnorm(x, birth.mean, birth.sd), min(xrange), max(xrange), add = TRUE,
lwd=3, col="blue")
I=c(1,5,25,50,75,95,99)
x1=qnorm(I/100, birth.mean,birth.sd)
y1=area*dnorm(x1,birth.mean, birth.sd)
label=c("1%", "5%", "25%", "50%","75%","95%","99%")
bit=5000
for(i in 1:7){
lines(c(x1[i],x1[i]), c(0,y1[i]),lwd=2,lty=2)
text(x1[i],-bit,label[i],cex=0.6)
}
lines(c(xw,xw),yrange,col="red",lwd=2)
# (d) Low birth weight
plot(xrange, yrange, type = "n", xlab = "Birthweight (gms)", ylab = "",
bty="n",axes=F,,main="(d) Low birth weight" )
axis(1)
curve(area*dnorm(x, birth.mean, birth.sd), min(xrange), max(xrange), add = TRUE,
lwd=3, col="blue")
x1=seq(1500,2500,10)
y1=area*dnorm(x1,birth.mean, birth.sd)
polygon(c(x1,x1[101:1]),c(rep(0,101), y1[101:1]),col="lightblue")
lines(c(xw,xw),yrange,col="red",lwd=2)
x1=seq(1500,xw,10)
nx=length(x1)
y1=area*dnorm(x1,birth.mean, birth.sd)
polygon(c(x1,x1[nx:1]),c(rep(0,nx), y1[nx:1]),col="red",density=10)
text(2000,70000,"Proportion\n below\n 2500 gms\n = 1.7%",cex=0.75)
text(3400,70000,"Proportion\n below\n 2910 gms\n = 11%",cex=0.75)
```
Entonces, podemos pensar en la población como un conjunto de individuos que provee
la distirbución de probabilidad de una observación aleatoria. Esto será muy útil
cuando lleguemos al momento de hacer *inferencia estadística.* Cuyo objetivo es
poder hacer afirmaciones sobre las características como media, moda, o dispersión que
en general no sabemos de antemano.
En los casos donde no hay muestreo (análisis de crímenes en una ciudad, o estudios
censales) la diferencia entre población y muestra no existe. Sin embargo, la noción
de población es valiosa. Pero, ¿cómo definimos una población?
```{block, type='comentario'}
Existen tres tipos de población de la cual podemos extraer una muestra de forma aleatoria:
- Población *literal.* Cuando podemos identificar a un grupo de dónde extraer muestras aleatorias.
- Población *virtual.* Cuando tomamos observaciones del ambiente, por ejemplo tomar mediciones de la calidad de aire. Los datos generados por en este escenario se denominan **muestras observacionales.**
- Población *metafórica.* Cuando no hay grupo de individuos mas grande. Pero aún asi podemos pensar como si los datos provienen de un espacio imaginario de posibilidades. Los datos geberados en este escenario se denominan **muestras naturales.**
```
## Trabajando como **Juez** {-}
Ahora nos enfocaremos en interpretar resultados estadísticos. Para esto
consideraremos el contexto de **dos muestras**. Es una práctica común en el análisis
estadístico pues permite contrastar el efecto de un diseño o prueba. Ejemplos
de esto los vemos al medir la tasa de captura en diseño de páginas *web*, pruebas
de una nueva medicina, o simple contraste entre dos poblaciones con
características distintas.
Una *respuesta* a una pregunta de interés
viene acomapañada de una *medida de incertidumbre*, la cual se basa en una
*modelo de probabilidad*. Hay veces en el que el modelo de probabilidad está
bien justificado, pues hay un *mecanismo aleatorio* detrás ---pensemos en el
lanzamiento de una moneda. En otras ocasiones el modelo de probabilidad es un
*artefacto* matemático que asemeja la realidad y permite aplicar
*modelos estadísticos*.
Para entender y comunicar las conclusiones de un modelo hay que estar
conscientes del mecanismo aleatorio que se utilizó, por ejemplo, en la selección
de unidades muestrales, o del grupo al que pertenecen.
Hay dos formas de hacer inferencia. La **inferencia causal** y la **inferencia a
poblaciones**. Saber los mecanismos que generaron los datos nos permite saber
qué tipo de inferencia es más adecuada para el estudio en cuestión.
### Inferencia Causal {-}
En un **experimento aleatorizado** la investigadora asume el control de
asignación de cada unidad experimental a los distintos grupos de estudio por
medio de un mecanismo aleatorio, por ejemplo, una moneda.
En un **estudio observacional** la asignación a los grupos se encuentra fuera del
control de la investigadora.
Es natural cuestionar si por medio de análisis estadísticos podemos concluir
relaciones causales. La respuesta es:
```{block, type='comentario'}
Las relaciones de causa y efecto se pueden inferir sólo si se utiliza un estudio aleatorizado,
pero no por medio de estudios observacionales.
```
El componente aleatorio asegura que las unidades observacionales con diferentes
características se mezclen, y cualquier evidencia de dicha relación se muestra
en el estudio. Aún asi, no hay certeza absoluta de la presencia de la
relación causal. Dicha incertidumbre es la que usualemente se pretende inculuir
en el modelo a través de técnicas estadísticas.
En un estudio observacional es imposible concluir una relación causal por medio
de un análisis estadístico. La analista no puede asegurar la ausencia de algún
factor de confusión (*confounding variable*) que sea responsable de distorsionar
las conclusiones.
```{block, type='comentario'}
Un factor de confusión está asociado tanto a la pertenencia de un grupo de estudio
como al resultado del estudio mismo. La presencia de un factor de confusión no permite
relacionar de manera directa la consecuencia con la pertenencia al grupo.
```
#### El valor de estudios observacionales {-}
Incluso aunque no podamos establecer relaciones causa-efecto, los estudios observacionales
poseen valor en un estudio formal. Las ventajas se pueden resumir en:
1. El objetivo del estudio. A veces establecer relaciones de causa-efecto no es el objetivo
2. Establecer la relacion causa-efecto se puede hacer por medio de otras rutas.
3. Datos observacionales pueden sugerir nuevas direcciones de investigación a través de *evidencia*.
#### Ejemplo: Policías y tráfico {-}
Supongamos que nos preguntan en cuánto reduce un policía el tráfico en
un crucero grande de la ciudad. La cultura popular
ha establecido que los policías en cruceros hacen más tráfico porque
no saben mover los semáforos.
Nosotros decidimos buscar datos para entender esto. Escogemos
entonces un grupo de cruceros problemáticos, registramos el tráfico
cuando visitamos, y si había un policía o no.
Después de este esfuerzo, obtenemos los siguientes datos:
```{r, message = FALSE, echo=FALSE}
library(tidyverse)
source("R/funciones_auxiliares.R")
n <- 5000
set.seed(881)
trafico_tbl <- tibble(x_inicial = rexp(n, 1 / 5), persistencia = runif(n)) %>%
mutate(dia_sem = sample(1:3, n, replace = T)) %>%
mutate(z = x_inicial + 10*persistencia - 15) %>%
mutate(policia = rbinom(n, 1, prob = exp(z) / (1 + exp(z)))) %>%
mutate(e = 0.2*rnorm(n)) %>%
mutate(y_0 = 1* x_inicial*exp(1.2 * persistencia - 1 + e)) %>%
mutate(y_1 = 1*x_inicial*exp(1.2 * persistencia - 1 - 0.6 * policia + e)) %>%
#mutate(y_0 = round(y_0, 3), y_1 = round(y_1)) %>%
mutate(tiempo_espera_min = round(ifelse(policia == 1, y_1, y_0), 2)) %>%
mutate(factor = x_inicial*exp(1.2 * persistencia - 1)) %>%
mutate(categoria = cut(factor, breaks = c(0, 4, 10, Inf),
labels = c("Fluido", "Típico", "Complicado"),
include.lowest = TRUE)) %>%
mutate(efecto = y_1 - y_0)
muestra_policias <- sample_n(trafico_tbl, 200) %>%
select(policia, tiempo_espera_min, categoria)
muestra <- muestra_policias %>% group_by(policia) %>% sample_n(5)
muestra %>% select(-categoria)
```
Lo que sabemos ahora es que la presencia de un policía es indicador
de tráfico alto. El análisis prosiguiría calculando medias y medidas de error
(escogimos una muestra aleatoria):
```{r, echo = FALSE, out.width ='70%', fig.align='center', fig.height =4, fig.asp =.45}
muestra %>% group_by(policia) %>%
summarise(media = mean(tiempo_espera_min), des_est = sd(tiempo_espera_min),
error = 2 * des_est / sqrt(n())) %>%
select(policia, media, error) %>%
mutate(policia = factor(policia)) %>%
mutate_if(is.numeric, round) %>%
ggplot(aes(x = policia, y = media, ymin = media - error, ymax = media + error)) +
geom_linerange() + coord_flip() + geom_point() +
ylab("Tiempo de espera promedio, minutos") +
labs(title = "Comparativo de tiempo de espera (mins)",
subtitle = "Intervalos de 95%" )
```
Si somos ingenuos, entonces podríamos concluir que los policías efectivamente
empeoran la situación cuando manipulan los semáforos, y confirmaríamos la
sabiduría popular.
Para juzgar este argumento desde el punto de vista causal, nos preguntamos primero:
- ¿Cuáles son los contrafactuales (los contrafactuales explican qué pasaría si hubiéramos
hecho otra cosa que la que efectivamente hicimos)
de las observaciones?
#### El estimador estándar {-}
A la comparación anterior ---la diferencia de medias de tratados y no tratados--- le llamamos usualmente el _estimador estándar_ del efecto causal. Muchas veces este es un estimador malo del efecto causal.
En nuestro ejemplo, para llegar a la conclusión errónea que confirma la sabiduría popular, hicimos un supuesto importante:
- En nuestra muestra, los casos con policía actúan como contrafactuales de los casos sin policía.
- Asi que asumimos que los casos con policía y sin policía son similares, excepto por la existencia o no de policía.
En nuestro ejemplo, quizá un analista más astuto nota que tienen
categorías históricas de qué tan complicado es cada crucero. Con esos datos obtiene:
```{r, echo = FALSE}
muestra
```
El analista argumenta entonces que los policías se enviaron principalmente a cruceros que
se consideran _complicados_ según datos históricos. Esto resta credibilidad a la
comparación que hicimos inicialmente:
- La comparación del estimador estándar no es de peras con peras: estamos comparando qué efecto tienen los
policías en cruceros difíciles con cruceros no difíciles donde no hay policía.
- La razón de esto es que el proceso generador de los datos incluye el hecho de que no
se envían policías a lugares donde no hay tráfico.
- ¿Cómo producir contrafactuales hacer la comparación correcta?
#### Experimentos tradicionales {-}
Idealmente, quisiéramos observar un mismo crucero en las dos condiciones: con y sin policías. Esto no es posible.
En un experimento "tradicional", como nos lo explicaron en la escuela, nos
aproximamos a esto preparando dos condiciones idénticas, y luego alteramos cada una de ellas
con nuestra intervención. Si el experimento está bien hecho, esto nos da observaciones
en pares, y cada quien tiene su contrafactual.
La idea del experimiento tradicional es _controlar_ todos los factores
que intervienen en los resultados, y sólo mover el tratamiento para producir
los contrafactuales. Más en general, esta estrategia consiste en hacer
_bloques_ de condiciones, donde las condiciones son prácticamente idénticas dentro e cada bloque. Comparamos entonces unidades tratadas y no tratadas
dentro de cada bloque.
Por ejemplo, si queremos saber si el tiempo de caída libre es diferente para un objeto
más pesado que otro, prepararíamos dos pesos con el mismo tamaño pero de peso distinto. Soltaríamos los dos al mismo tiempo y compararíamos el tiempo de caída de cada uno.
En nuestro caso, como es usual en problemas de negocio o sociales, hacer esto es considerablemente más difícil. No podemos "preparar" cruceros con condiciones idénticas. Sin embargo, podríamos intentar bloquear los cruceros
según información que tenemos acerca de ellos, para hacer más comparaciones e peras con peras.
#### Bloqueo {-}
Podemos acercanos en lo posible a este ideal de experimentación usando
información existente.
En lugar de hacer comparaciones directas entre unidades que recibieron
el tratamiento y las que no (que pueden ser diferentes en otros
aspectos, como vimos arriba),
podemos refinar nuestras comparaciones _bloquéandolas_ con variables
conocidas.
En el ejemplo de los policías, podemos hacer lo siguiente: dentro de
_cada categoría de cruceros_ (fluido, típico o complicado), tomaremos una muestra de cruceros, algunos con
policía y otros sin. Haremos comparaciones dentro de cada categoría.
Obtenemos un muestra con estas características (6 casos en cada categoría
de crucero, 3 con policía y 3 sin policía):
```{r, echo = FALSE}
muestra_bloqueada <- trafico_tbl %>% group_by(categoria, policia) %>%
sample_n(3) %>%
select(policia, tiempo_espera_min, categoria)
knitr::kable(muestra_bloqueada %>% group_by(categoria, policia) %>% tally())
```
Y ahora hacemos comparaciones dentro de cada bloque creado por categoría:
```{r, message = FALSE, echo = FALSE}
muestra_bloqueada %>%
group_by(categoria, policia) %>%
summarise(tiempo_espera = mean(tiempo_espera_min)) %>%
mutate_if(is.numeric, ~ round(.x, 1)) %>%
pivot_wider(names_from = policia,
values_from = tiempo_espera,
names_prefix = "policia =")
```
Y empezamos a ver otra imagen en estos datos: comparando tipos
e cruceros similares, los que tienen policía tienen tiempos de
espera ligeramente más cortos.
¿Hemos termniado? ¿Podemos concluir que el efecto de un policía
es beneficiosos pero considerablemente chico? ¿Qué problemas
puede haber con este análisis?
#### Variables desconocidas {-}
El problema con el análisis anterior es que controlamos por una
variable que conocemos, pero muchas otras variables pueden estar
ligadas con el proceso de selección de cruceros para enviar policías.
- Por ejemplo, envían o policías a cruceros _Típicos_ solo cuando
reportan mucho tráfico.
- No envían a un polícia a un crucero _Complicado_ si no presenta demasiado
tráfico.
- Existen otras variables desconocidas que los tomadores de decisiones
usan para enviar a los policías.
En este caso, por ejemplo, los expertos hipotéticos
nos señalan que hay algunos
cruceros que aunque problemáticos a veces, su tráfico se resuelve
rápidamente, mientras que otros tienen tráfico más persistente, y
prefieren enviar policías a los de tráfico persistente. La lista
de cruceros persistentes están en una hoja de excel que se comparte
de manera informal.
En resumen, no tenemos conocimiento detallado del **proceso generador
de datos** en cuanto a cómo se asignan los policías a los cruceros.
Igual que en la sección anterior, podemos cortar esta complejidad
usando **aleatorización**.
Nótese que los expertos no están haciendo nada malo: en su trabajo
están haciendo el mejor uso de los recursos que tienen. El problema
es que por esa misma razón no podemos saber el resultado de sus esfuerzos,
y si hay maneras de optimizar la asignación que hacen actualmente.
#### Aleatorizando el tratamiento {-}
Tomamos la decisión entonces de hacer un experimento que incluya
aletorización.
En un día
particular, escogeremos algunos cruceros.
Dicidimos usar solamente cruceros de la categoría _Complicada_ y
_Típica_, pues
esos son los más interesantes para hacer intervenciones.
Usaremos un poco de código para entener el detalle: en estos datos,
tenemos para cada caso los dos posibles resultados ipotéticos
$y_0$ y $y_1$ (con
policia y sin policia). En el experimento asignamos el
tratamiento al azar:
```{r}
muestra_exp <- trafico_tbl %>% filter(categoria != "Fluido") %>%
sample_n(200) %>%
# asignar tratamiento al azar, esta es nuestra intervención:
mutate(tratamiento_policia = rbernoulli(length(y_0), 0.5)) %>%
# observar resultado
mutate(tiempo_espera_exp = ifelse(tratamiento_policia ==1, y_1, y_0))
```
Nótese la diferencia si tomamos la asignación natural del tratamiento (policía o no):
```{r}
set.seed(134)
muestra_natural <- trafico_tbl %>% filter(categoria != "Fluido") %>%
sample_n(200) %>%
# usamos el tratamiento que se asignó
# policia indica si hubo o no policía en ese crucero
# observar resultado
mutate(tiempo_espera_obs = ifelse(policia ==1, y_1, y_0))
```
Resumimos nuestros resultados del experimento son:
```{r, echo = FALSE}
muestra_exp %>%
mutate(tratamiento_policia = as.numeric(tratamiento_policia)) %>%
group_by(categoria, tratamiento_policia) %>%
summarise(tiempo_espera = mean(tiempo_espera_exp)) %>%
pivot_wider(names_from = tratamiento_policia,
values_from = tiempo_espera,
names_prefix = "policia=")
```
Sin embargo, la muestra natural da:
```{r, echo = FALSE}
muestra_natural %>% group_by(categoria, policia) %>%
summarise(tiempo_espera = mean(tiempo_espera_obs)) %>%
pivot_wider(names_from = policia,
values_from = tiempo_espera,
names_prefix = "policia=")
```
**¿Cuál de los dos análisis da la respuesta correcta a la pregunta:
ayudan o no los policías a reducir el tráfico en los cruceros
problemáticos?** El experimento establece que un policía en promedio
reduce a la mitad el tiempo de espera en un crucero complicado.
### Inferencia a poblaciones {-}
La situación es bastante clara. Inferir características de una poblacion **sólo** se puede realizar por
medio de muestreo aleatorio, no de otra forma.
Seleccionar de manera aleatoria significa que cualquier conjunto de tamaño $N$ que escojamos tiene
la misma probabilidad de ser escogido que cualquier otro conjunto del mismo tamaño.
### Resumen: selección de unidades y tratamiento {-}
Vimos dos tipos de inferencia que requieren distintos diseños de estudio,
en particular debemos considerar el mecanismo de aleatorización para
entender las inferencias que podemos hacer: causal o a poblaciones.
El punto crucial para entender las medidas de incertidumbre estadística es
visualizar de manera hipotética, replicaciones del estudio y las condiciones
que llevaron a la selección de la muestra. Esto es, entender el proceso
generador de datos e imaginar replicarlo.
![Inferencia estadística de acuerdo al tipo del diseño [@ramsey]](images/03_inferencia-estudio.png)
* El cuadro en la esquina superior izquierda es donde el análisis es más simple y los
resultados son más fáciles de interpretar.
* Es posible hacer análisis fuera de este cuadro, pero el proceso es más
complicado, requieren más supuestos, conocimiento del dominio y habilidades
de análisis. En general resultan conclusiones menos sólidas. Muchas veces no
nos queda otra opción más que trabajar fuera del cuadro ideal.
```{block, type='ejercicio'}
Ubica los siguientes tipos de análisis:
- Pruebas clínicas para medicinas
- Analizar cómo afecta tener seguro médico a los ingresos, usando datos del ENIGH.
- Estimación de retorno sobre inversión en modelos de marketing mix.
```
#### Asignación natural del tratamiento {-}
- Cuando consideramos un sistema donde se "asignan" tratamientos,
generalmente los tratamientos se asignan bajo un criterio de
optimización o conveniencia.
- La cara buena de este hecho es que de alguna forma los resultados
están intentando optimizarse, y la gente está haciendo su trabajo.
- La cara mala de este hecho es que no podemos evaluar de manera simple la
efectividad de los tratamientos. Y esto hace difícil **optimizar** de forma
cuantificable los procesos, o **entender** qué funciona y qué no.
<!-- ### Experimentos y datos observacionales {-} -->