diff --git a/02-tipos-de-estudio.Rmd b/02-tipos-de-estudio.Rmd index 11244c1..508433f 100644 --- a/02-tipos-de-estudio.Rmd +++ b/02-tipos-de-estudio.Rmd @@ -36,6 +36,21 @@ ggplot(bind_rows(pre, post), aes(x = y, color = group)) + xlab("CPR error") + labs(color = "") ``` + + + + + + + + + + + + + + +
> Far better an approximate answer to the right question, which is often vague, than an exact answer to the wrong question, which can always be made precise. @@ -84,7 +99,8 @@ ggplot(puntos, aes(x = x, y = y)) + Ejemplos: Alguien nos pregunta cuáles son las tiendas que mas venden de una cadena. Podríamos consultar bases de datos, hacer extracciones, definir -periodos, etc. y dar una respuesta que probablemente es poco útil. Nos damos +periodos, etc. y reportar el promedio de ventas en el último mes, esta respuesta +probablemente es poco útil. Nos damos cuenta, por ejemplo, porque la peor tienda es una que abrió hace relativamente poco, y la mejor es una de las tiendas más grandes que está en una zona de tráfico de alto costo. Una pregunta más interesante @@ -118,7 +134,8 @@ decimos que puede ocurrir con 10% de probabilidad ocurre efectivamente 1 de cada 10 veces, si decimos 20% entonces ocurre 2 de 20, etc. Veremos que para muestras dadas naturalmente, a veces es muy difiícil entender a -fondo el proceso generación de la muestra. +fondo el proceso que generó la muestra y por tanto no tenemos las garantías de +eficiencia y calibración. ### Ejemplo: Prevalencia de anemia {-} @@ -176,6 +193,10 @@ 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). +El elemento clave, es la aleatorización en la selección de la muestra, la idea +es distribuir los efecros desconcidos o no controlables que pueden introducir +sesgos o variabilidad no conocida en los resultados. + ## Pero si no podemos hacer muestreo aleatorio? {-} En el caso de prevalencia de anemia, discutiendo con médicos e investigadores @@ -206,8 +227,6 @@ 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? - ### Ejemplo: Policías y tráfico {-} @@ -244,9 +263,9 @@ trafico_tbl <- tibble(x_inicial = rexp(n, 1 / 5), include.lowest = TRUE)) |> mutate(efecto = y_1 - y_0) muestra_policias <- sample_n(trafico_tbl, 200) |> - select(policia, tiempo_espera_min, categoria) + dplyr::select(policia, tiempo_espera_min, categoria) muestra <- muestra_policias |> group_by(policia) |> sample_n(5) -muestra |> select(-categoria) +muestra |> dplyr::select(-categoria) ``` Lo que sabemos ahora es que la presencia de un policía es indicador @@ -309,7 +328,7 @@ se consideran _Complicados_ según datos históricos. Esto resta credibilidad a 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. +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 para hacer la comparación correcta? @@ -356,7 +375,7 @@ 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) + dplyr::select(policia, tiempo_espera_min, categoria) knitr::kable(muestra_bloqueada |> count(categoria, policia)) ``` @@ -488,7 +507,9 @@ fundamental para entender las inferencias que podemos hacer en distintos escenar ![Inferencia estadística de acuerdo al tipo del diseño (@ramsey).](images/03_inferencia-estudio.png) * El cuadro arriba a la izquierda es donde el análisis es más simple y los -resultados son más fáciles de interpretar. +resultados son más fáciles de interpretar. En este escenario don de la aleatorización +es tanto en unidades como en grupos no hacen falta supuestos adicionales para +tener las garantías de métodos de inferencia. * 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 diff --git a/04-distribucion-muestreo.Rmd b/04-distribucion-muestreo.Rmd new file mode 100644 index 0000000..906144b --- /dev/null +++ b/04-distribucion-muestreo.Rmd @@ -0,0 +1,756 @@ +# Estimación y distribución de muestreo + +```{r setup, include=FALSE, message=FALSE} +library(tidyverse) +source("R/funciones_auxiliares.R") +knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning=FALSE, fig.align = 'center',fig.width = 5, fig.height=3) +comma <- function(x) format(x, digits = 2, big.mark = ",") +theme_set(theme_minimal()) +``` + +En esta sección discutiremos cuál el objetivo general del proceso de estimación. +y cómo entender y manejar la variabilidad que se produce cuando aleatorizamos +la selección de las muestras que utilizamos para hacer análisis. + +A diferencia de las pruebas de permutación, donde evaluábamos como cambiaría una +estadísitica si un tratamiento o grupo se hubiera asignado de forma distinta, en la +siguiente sección nos preguntamos como varía una estadística entre muestras. Por ejemplo, +pasaremos de preguntar si una vacuna reduce el riesgo de una enfermedad a evaluar +en que magnitud se reduce el riesgo de contraer la enfermedad. + + + +## Ejemplo: precios de casas {-} + +Supongamos que queremos conocer el valor total de las casas +que se vendieron recientemente en una zona +particular. Supondremos que tenemos un listado de las casas que se han +vendido recientemente, pero en ese listado no se encuentra el precio de venta. +Decidimos entonces tomar una muestra aleatoria de 100 de esas casas. Para esas +casas hacemos trabajo de campo para averiguar el precio de venta. + +```{r} +marco_casas <- read_csv("data/casas.csv") +set.seed(841) +muestra_casas <- sample_n(marco_casas, 100) |> + select(id, nombre_zona, area_habitable_sup_m2, precio_miles) +sprintf("Hay %0.0f casas en total, tomamos muestra de %0.0f", + nrow(marco_casas), nrow(muestra_casas)) +head(muestra_casas) +``` + +Como tomamos una muestra aleatoria, intentamos estimar el valor +total de las casas que se vendieron expandiendo el total muestral, es decir +nuestro estimador $\hat{t} = t(X_1,\ldots X_{100})$ del total +poblacional $t$ es + +$$\hat{t} = \frac{N}{n} \sum_{i=1}^{100} X_i = N\bar{x}$$ + +Esta función implementa el estimador: + +```{r} +n <- nrow(muestra_casas) # tamaño muestra +N <- nrow(marco_casas) # tamaño población +estimar_total <- function(muestra_casas, N){ + total_muestral <- sum(muestra_casas$precio_miles) + n <- nrow(muestra_casas) + # cada unidad de la muestra representa a N/n + f_exp <- N / n + # estimador total es la expansión del total muestral + estimador_total <- f_exp * total_muestral + res <- tibble(total_muestra = total_muestral, + factor_exp = f_exp, + est_total_millones = estimador_total / 1000) + res +} +estimar_total(muestra_casas, N) |> + mutate(across(where(is.numeric), \(x) round(x, 2))) +``` + +Sin embargo, si hubiéramos obtenido otra muestra, hubiéramos obtenido otra +estimación diferente. Por ejemplo: + +```{r} +estimar_total(sample_n(marco_casas, 100), N) |> + mutate(across(where(is.numeric), round, 2)) +``` + +El valor poblacional que buscamos estimar (nótese que en la práctica este no lo conocemos) +es: + +```{r} +# multiplicar por 1000 para que sea en millones de dólares +total_pob <- sum(marco_casas |> pull(precio_miles)) / 1000 +total_pob +``` + +Así que: + +- Para algunas muestras esta estadística puede estar muy cercana al valor poblacional, +pero para otras puede estar más lejana. +- Para entender qué tan buena es una estimación +particular, entonces, tenemos que entender *cuánta variabilidad hay de muestra a muestra* +debida a la aleatorización. Esto depende del diseño de la muestra y +de la población de precios de casas (que no conocemos). + +## Distribución de muestreo {-} + +La distribución de muestreo de una estadística enumera los posibles resultados +que puede tomar esa estadística sobre todas las muestras posibles. Este es el concepto +básico para poder entender qué tan bien o mal estima un parámetro poblacional dado. + +En nuestro ejemplo anterior de precio de casas, no podemos calcular todas las posibles +estimaciones bajo todas las posibles muestras, pero podemos aproximar +repitiendo una gran cantidad de veces el proceso de muestreo, como hicimos +al aproximar la distribución de permutaciones de estadísticas de prueba de las +secciones anteriores. + +Empezamos repitiendo 10 veces y examinamos cómo varía nuestra estadística: + +```{r} +replicar_muestreo <- function(marco_casas, m = 500, n){ + # n es el tamaño de muestra que se saca de marco_casas + # m es el número de veces que repetimos el muestro de tamaño n + resultados <- map_df(1:m, + function(id) { + sample_n(marco_casas, n) |> + estimar_total(N) + }, .id = "id_muestra") +} +replicar_muestreo(marco_casas, m = 10, n = 100) |> + mutate(across(where(is.numeric), round, 1)) |> + formatear_tabla() +``` + +Como vemos, hay variación considerable en nuestro estimador del total, pero +la estimación que haríamos con cualquiera de estas muestras no es muy mala. Ahora +examinamos un número más grande de simulaciones: + +```{r, cache = TRUE} +replicaciones_1 <- replicar_muestreo(marco_casas, m = 1500, n = 100) +``` + +Y el siguiente histograma nos dice qué podemos esperar de la variación de +nuestras estimaciones, y donde es más probable que una estimación particular caiga: + +```{r} +graf_1 <- ggplot(replicaciones_1, aes(x = est_total_millones)) + + geom_histogram() + + geom_vline(xintercept = total_pob, colour = "red") + + xlab("Millones de dólares") + + scale_x_continuous(breaks = seq(180, 240, 10), limits = c(180, 240)) +graf_1 +``` + +Con muy alta probabilidad el error no será de más de unos 30 millones de dólares +(o no más de 20% del valor poblacional). + + +```{block, type='mathblock'} +**Definición** Sea $X_1, X_2, \ldots X_n$ una muestra, +y $T = t(X_1, X_2, \ldots, X_n)$ una estadística. + +La **distribución de muestreo** de $T$ es la función de distribución de $T$. Esta +distribución es sobre todas las posibles muestras que se pueden obtener. + +Cuando usamos $T$ para estimar algún parámetro poblacional $\theta$, decimos +informalmente que el estimador es **preciso** si su distribución de muestreo está +muy concentrada alrededor del valor $\theta$ que queremos estimar. +``` + +Si la distribución de muestreo está concentrada en un conjunto muy grande +o muy disperso, quiere decir que con alta probabilidad cuando obtengamos +nuestra muestra y calculemos nuestra estimación, el resultado estará lejano +del valor poblacional que nos interesa estimar. + +Veamos qué pasa cuando hacemos la muestra más grande en nuestro ejemplo: + +```{r, cache = TRUE} +replicaciones_2 <- replicar_muestreo(marco_casas, m = 1500, n = 250) +``` + +Graficamos las dos distribuciones de muestreo juntas, y vemos cómo +con mayor muestra obtenemos un estimador más preciso, y sin considerar el costo, +preferimos el estimador **mejor concentrado alrededor del valor que buscamos estimar**. + +```{r, fig.width = 8, fig.height = 4} +library(patchwork) +graf_2 <- ggplot(replicaciones_2, aes(x = est_total_millones)) + + geom_histogram() + + geom_vline(xintercept = total_pob, colour = "red") + + xlab("Millones de dólares") + + scale_x_continuous(breaks = seq(180, 240, 10), limits = c(180, 240)) +graf_1 + graf_2 +``` + +```{block, type = 'comentario'} +**Observación**: a veces este concepto se confunde la distribución +poblacional de las $X_i$. Esto es muy diferente. + +``` + +Por ejemplo, en nuestro caso, el histograma de la distribución de valores poblacionales es + +```{r, fig.width = 5, fig.height=3} +ggplot(marco_casas, aes(x = precio_miles)) + geom_histogram() +``` +que en general no tiene ver mucho en escala o forma con la distribución de muestreo +de nuestro estimador del total. + +## Más ejemplos {-} + +Podemos también considerar muestrear de poblaciones sintéticas o modelos +probabilísticos que usamos para modelar poblaciones reales. + +Por ejemplo, supongamos que tomamos una muestra de tamaño 15 de la distribución +uniforme en $[0,1]$. Es decir, cada $X_i$ es un valor uniformemente distribuido +en $[0,1]$, y las $X_i$ se extraen independientemente unas de otras. Consideramos +dos estadísticas de interés: + +1. La media muestral $T_1(X) = \frac{1}{15}\sum_{i = 1}^{15} X_i$ +2. El cuantil 0.75 de la muestra $T_2(X) = q_{0.75}(X)$ + +```{block2, type='ejercicio'} +¿Cómo crees que se vean las distribuciones muestrales de estas estadísticas? + ¿Alrededor de qué valores crees que concentren? ¿Crees que tendrán mucha o poca +dispersión? ¿Qué forma crees que tengan? +``` + +Para el primer caso hacemos: + +```{r} +# simular +replicar_muestreo_unif <- function(est = mean, m, n = 15){ + valores_est <- map_dbl(1:m, ~ est(runif(n))) + tibble(id_muestra = 1:m, estimacion = valores_est) +} +sim_estimador_1 <- replicar_muestreo_unif(mean, 4000, 15) +# graficar aprox de distribución de muestreo +ggplot(sim_estimador_1, aes(x = estimacion)) + + geom_histogram(bins = 40) + + xlim(c(0, 1)) +``` + +```{r} +# simular para el máximo +cuantil_75 <- function(x) quantile(x, 0.75) +sim_estimador_2 <- replicar_muestreo_unif(cuantil_75, 4000, 15) +# graficar distribución de muestreo +ggplot(sim_estimador_2, aes(x = estimacion)) + + geom_histogram(breaks = seq(0, 1, 0.02)) + + xlim(c(0, 1)) +``` + + +```{block2, type='ejercicio'} +Supón que tenemos una muestra de 30 observaciones de una distribución +uniforme $[0,b]$. + +- ¿Qué tan buen estimador de $b/2$ es la media muestral? ¿Cómo lo cuantificarías? +- ¿Qué tan buen estimador del cuantil 0.8 de la distribución uniforme es +el cuantil 0.8 muestral? ¿Qué desventajas notas en este estimador? + +``` + +## El error estándar {-} + +Una primera medida útil de la dispersión de la distribución de muestreo +es su desviación estándar: la razón específica tiene qué ver +con un resultado importante, el teorema central del límite, que veremos +más adelante. En este caso particular, a esta desviación estándar +se le llama error estándar: + +```{block, type='mathblock'} +**Definición** A la desviación estándar de una estadística $T$ le llamamos +su **error estándar**, y la denotamos por $\text{ee}(T)$. A cualquier estimador +de este error estándar lo denotamos como $\hat{\text{ee}}(T)$. + +``` + +Este error estándar mide qué tanto varía el estimador $T$ de muestra a muestra. + +```{block, 'comentario'} +**Observación**: es importante no confundir el error estándar con +la desviación estándar de una muestra (o de la población). +``` + + +En nuestro ejemplo +de las uniformes, la desviación estándar de las muestras varía como: + +```{r} +map_dbl(1:10000, ~ sd(runif(15))) |> quantile() |> round(2) +``` + +Mientras que el error estándar de la media es aproximadamente + +```{r} +map_dbl(1:10000, ~ mean(runif(15))) |> sd() +``` + +y el error estándar del máximo es aproximadamente + +```{r} +map_dbl(1:10000, ~ max(runif(15))) |> sd() +``` + +```{block2, type='ejercicio'} +Como ejercicio para contrastar estos conceptos, +puedes considerar: ¿Qué pasa con la desviación estándar de una muestra muy +grande de uniformes? ¿Qué pasa con el error estándar de la media muestral de una muestra muy grande de uniformes? +``` + + + +### Ejemplo: valor de casas {-} + +Consideramos el error estándar del estimador del total del inventario vendido, usando +una muestra de 250 con el estimador del total que describimos arriba. Como aproximamos con +simulación la distribución de muestreo, podemos hacer: + +```{r} +ee_2 <- replicaciones_2 |> pull(est_total_millones) |> sd() +round(ee_2, 1) +``` + +que está en millones de pesos y cuantifica la dispersión de la distribución de +muestreo del estimador del total. + +Para tamaño de muestra 100, obtenemos más dispersión: + +```{r} +ee_1 <- replicaciones_1 |> pull(est_total_millones) |> sd() +round(ee_1, 1) +``` + +Nótese que esto es muy diferente, por ejemplo, a la desviación estándar +poblacional o de una muestra. Estas dos cantidades miden la variabilidad del +estimador del total. + +## Calculando la distribución de muestreo {-} + +En los ejemplos anteriores usamos simulación para obtener aproximaciones +de la distribución de muestreo de algunos estimadores. También +es posible: + +- Hacer cálculos exactos a partir de modelos +probabilísticos. +- Hacer aproximaciones asintóticas para muestras grandes (de las cuales +la más importante es la que da el teorema central del límite). + +En los ejemplos de arriba, cuando muestreamos de la poblaciones, +extrajimos las muestras de manera aproximadamente independiente. Cada +observación $X_i$ tiene la misma distribución y las $X_i$'s son +independientes. Este tipo de diseños aleatorizados es de los más +simples, y se llama **muestreo aleatorio simple**. + +En general, en esta parte haremos siempre este supuesto: Una **muestra** +es iid (independiente e idénticamente distribuida) si es +es un conjunto de observaciones $X_1,X_2, \ldots X_n$ independientes, +y cada una con la misma distribución. + +En términos de poblaciones, esto lo logramos obteniendo cada observación +de manera aleatoria con el mismo procedimiento. En términos de modelos +probabilísticos, cada $X_i$ se extrae de la misma distribución fija $F$ +(que pensamos como la "población") de manera independiente. Esto lo denotamos +por $X_i \overset{iid}{\sim} F.$ + +### Ejemplo {-} + +Si $X_1, X_2, \ldots X_n$ es una muestra de uniformes independientes en $[0,1]$, ¿cómo +calcularíamos la distribución de muestreo del máximo muestra $T_2 = \max$? En este +caso, es fácil calcular su función de distribución acumulada de manera exacta: + +$$F_{\max}(x) = P(\max\{X_1,X_2,\ldots X_n\} \leq x)$$ +El máximo es menor o igual a $x$ si y sólo si todas las $X_i$ son menores +o iguales a $x$, así que +$$F_{\max} (x) = P(X_1\leq x, X_2\leq x, \cdots, X_n\leq x)$$ +como las $X_i$'s son independientes entonces +$$F_{\max}(x) = P(X_1\leq x)P(X_2\leq x)\cdots P(X_n\leq x) = x^n$$ +para $x\in [0,1]$, pues para cada $X_i$ tenemos $P(X_i\leq x) = x$. +Así que no es necesario usar simulación para conocer esta distribución de muestreo. +Derivando esta distribución acumulada obtenemos su densidad, que es + +$$f(x) = nx^{n-1}$$ + +para $x\in [0,1]$, y es cero en otro caso. + +Si comparamos con nuestra simulación: + +```{r} +teorica <- tibble(x = seq(0, 1 ,0.001)) |> + mutate(f_dens = 15 * x^14) +sim_estimador_3 <- replicar_muestreo_unif(max, 4000, 15) +ggplot(sim_estimador_3) + + geom_histogram(aes(x = estimacion), breaks = seq(0, 1, 0.02)) + + xlim(c(0.5, 1)) + + # el histograma es de ancho 0.02 y el número de simulaciones 4000 + geom_line(data = teorica, aes(x = x, y = (4000 * 0.02) * f_dens), + colour = "red", linewidth = 1.3) +``` + +Y vemos que con la simulación obtuvimos una buena aproximación + + +**Nota**: ¿cómo se relaciona un histograma con la función de densidad +que genera los datos? Supón que $f(x)$ es una función de densidad, y +obtenemos un número $n$ de simulaciones independientes. Si escogemos +un histograma de ancho $\Delta$, ¿cuántas observaciones esperamos +que caigan en un intervalo $I = [a - \Delta/2, a + \Delta/2]$?. La probabilidad +de que una observación caiga en $I$ es igual a + +$$P(X\in I) = \int_I f(x)\,dx = \int_{a - \Delta/2}^{a + \Delta/2} f(x)\,dx \approx f(a) \text{long}(I) = f(a) \Delta$$ +para $\Delta$ chica. Si nuestra muestra es de tamaño $n$, el número esperado +de observaciones que caen en $I$ es entonces $nf(a)\Delta$. Eso explica +el ajuste que hicimos en la gráfica de arriba. Otra manera de hacer es +ajustando el histograma: si en un intervalo el histograma alcanza el valor $y$, +$$f(a) = \frac{y}{n\Delta}$$ + +```{r} +teorica <- tibble(x = seq(0, 1 ,0.001)) |> + mutate(f_dens = 15*x^{14}) +ggplot(sim_estimador_3) + + geom_histogram(aes(x = estimacion, y = after_stat(density)), breaks = seq(0, 1, 0.02)) + + xlim(c(0.5, 1)) + + # el histograma es de ancho 0.02 y el número de simulaciones 4000 + geom_line(data = teorica, aes(x = x, y = f_dens), + colour = "red", size = 1.3) +``` + +### Ejemplo {-} + +Supongamos que las $X_i$'s son independientes y exponenciales con tasa $\lambda > 0$. +¿Cuál es la distribución de muestreo de la suma $S = X_1 + \cdots + X_n$? Sabemos que la suma +de exponenciales independientes es una distribución gamma con parámetros $(n, \lambda)$, +y esta es la distribución de muestreo de nuestra estadística $S$ bajo las hipótesis +que hicimos. + +Podemos checar este resultado con simulación, por ejemplo para una +muestra de tamaño $n=15$ con $\lambda = 1$: + +```{r} +replicar_muestreo_exp <- function(est = mean, m, n = 150, lambda = 1){ + valores_est <- map_dbl(1:m, ~ est(rexp(n, lambda))) + tibble(id_muestra = 1:m, estimacion = valores_est) +} +sim_estimador_1 <- replicar_muestreo_exp(sum, 4000, n = 150) +teorica <- tibble(x = seq(0, 35, 0.001)) |> + mutate(f_dens = dgamma(x, shape = 150, rate = 1)) +# graficar aprox de distribución de muestreo +ggplot(sim_estimador_1) + + geom_histogram(aes(x = estimacion, y = ..density..), bins = 35) + + geom_line(data = teorica, aes(x = x, y = f_dens), colour = "red", size = 1.2) +``` + +## Teorema central del límite {-} + +Si consideramos los ejemplos de arriba donde tratamos con estimadores +basados en una suma, total o una media ---y en menor medida cuantiles muestrales---, +vimos que las distribución de +muestreo de las estadísticas que usamos tienden a tener una forma común. +Estas son manifestaciones de una regularidad estadística importante que +se conoce como el **teorema central del límite**: las distribuciones de muestreo +de sumas y promedios son aproximadamente normales cuando el tamaño de muestra +es suficientemente grande. + +```{block2, type="mathblock"} +**Teorema central del límite** + + Si $X_1,X_2, \ldots, X_n$ son independientes e idénticamente distribuidas con +media $\mu$ y desviación estándar $\sigma$ finitas. + +Si el tamaño de muestra $n$ es grande, entonces la distribución de muestreo de la media + +$$\bar{X} = \frac{X_1 + X_2 +\cdots + X_n}{n}$$ + + es aproximadamente normal con media $\mu$ y desviación estándar $\sigma/\sqrt{n}$, +que escribimos como + +$$\bar{X} \xrightarrow{} \mathsf{N}\left( \mu, \frac{\sigma}{\sqrt{n}} \right)$$ + +Adicionalmente, la distribución de la +media estandarizada converge a una distribución normal +estándar cuando $n$ es grande: +$$\sqrt{n} \, \left( \frac{\bar{X}-\mu}{\sigma} \right) \xrightarrow{} \mathsf{N}(0, 1)$$ + +``` + +- El error estándar de $\bar{X}$ es +$\text{ee}(\bar{X}) = \frac{\sigma}{\sqrt{n}}$. Si tenemos una muestra, podemos +estimar $\sigma$ con de la siguiente forma: +$$\hat{\sigma} =\sqrt{\frac{1}{n}\sum_{i=1}^n (X_i - \bar{X})^2}$$ +o el más común (que explicaremos más adelante) +$$\hat{s} = \sqrt{\frac{1}{n-1}\sum_{i=1}^n (X_i - \bar{X})^2}$$ + +- Este hecho junto con el teorema del límite central nos dice cuál es la dispersión, +y cómo se distribuyen las posibles desviaciones de la media muestral alrededor +de la verdadera media poblacional. + +- ¿Qué tan grande debe ser $n$. Depende de cómo es la población. Cuando la población +tiene una distribución muy sesgada, por ejemplo, $n$ típicamente +necesita ser más grande que cuando la población es simétrica si queremos +obtener una aproximación "buena". + +- En algunos textos se afirma que $n\geq 30$ es suficiente para que la +aproximación del Teorema central del límite (TCL) sea buena siempre y cuando +la distribución poblacional no sea muy sesgada. Esta regla es más o menos arbitraria +y es mejor no confiarse, pues fácilmente puede fallar. En la práctica es importante +checar este supuesto, por ejemplo usando remuestreo (que veremos más adelante) + +```{block, type='ejercicio'} +Revisa los ejemplos que hemos visto hasta ahora (precios de casas, simulaciones +de uniformes y exponenciales según las distintas estadísticas que consideramos). +¿Qué distribuciones de muestreo +parecen tener una distribución normal? ¿Cómo juzgamos si estas distribuciones están +cerca o lejos de una distribución normal? +``` + +## Normalidad y gráficas de cuantiles normales {-} + +Para checar si una distribución de datos dada es similar a la normal, la herramienta +mas común en estádística es la gráfica de cuantiles teóricos, que es una generalización +de la gráfica de cuantiles que vimos anteriormente. + +En primer lugar, definimos la función de cuantiles de una distribución teórica, +que es análoga a la que definimos para conjuntos de datos: + +Supongamos que tenemos una distribución acumulada teórica $\Phi$. Podemos +definir el cuantil-$f$ $q(f)$ de $\Phi$ como +el valor $q(f)$ tal que +$$q(f) = \text{argmin}\{x \,| \, \Phi(x)\geq f \}$$ + +En el caso de que $\Phi$ tiene densidad $\phi$, y su soporte es un intervalo (que puede +ser de longitud infinita), entonces podemos también escribir $q(f)$ +como el valor único donde acumulamos $f$ de la probabilidad + +$$\int_{-\infty}^{q(f)} \phi(x)\,dx= f$$ +Por ejemplo, para una densidad normal, abajo mostramos los cuantiles $f=0.5$ (mediana) +y $f=0.95$ + +```{r} +densidad_tbl <- tibble(x = seq(0, 10, 0.01)) |> + mutate(densidad = dnorm(x, 5, 1)) +``` + +```{r, fig.width=7, fig.height=3} +# qnorm es la función de cuantiles de una normal +cuantil_50 <- qnorm(0.50, 5, 1) +cuantil_90 <- qnorm(0.95, 5, 1) +# graficamos +densidad_tbl <- densidad_tbl |> + mutate(menor_50 = x >= cuantil_50) |> + mutate(menor_90 = x >= cuantil_90) +g_normal_50 <- ggplot(densidad_tbl, aes(y = densidad)) + + ylab('f(x)') + + geom_area(aes(x = x, fill = menor_50)) + + geom_line(aes(x = x), alpha = 0.1) + + geom_vline(xintercept = cuantil_50) + theme(legend.position = "none") + + annotate("text", 4.3, 0.2, label = "50%") + + labs(subtitle = paste0("q(0.5)=", round(cuantil_50,1))) +g_normal_90 <- ggplot(densidad_tbl, aes(y = densidad)) + + ylab('f(x)') + + geom_area(aes(x = x, fill = menor_90)) + + geom_line(aes(x = x), alpha = 0.1) + + geom_vline(xintercept = cuantil_90) + theme(legend.position = "none") + + annotate("text", 5.0, 0.2, label = "95%") + + labs(subtitle = paste0("q(0.95)=", round(cuantil_90,1))) +g_normal_50 + g_normal_90 +``` + +Como todas las distribuciones normales tienen la misma forma, y para obtener una de otra +solo basta reescalar y desplazar, para calcular los cuantiles de una +variable con distribución normal $\mathsf{N}(\mu, \sigma)$ +sólo tenemos que saber los cuantiles de la distribución normal estándar $\mathsf{N}(0,1)$ y escalarlos +apropiadamente por su media y desviación estándar + +$$q(f, \mu, \sigma) = \mu + \sigma q(f, 0, 1)$$ +Puedes demostrar esto sin mucha dificultad empezando con $P(X\leq q) = f$ y estandarizando: + +$$P(X\leq q(f, \mu, \sigma)) = f \implies P\left (Z\leq \frac{q(f,\mu,\sigma) - \mu}{\sigma}\right)=f$$ +y esto implica que +$$q(f, 0, 1) = \frac{q(f,\mu,\sigma) - \mu}{\sigma} \implies q(f, \mu, \sigma) = \mu + \sigma q(f, 0, 1)$$ + +De modo que si graficáramos los cuantiles de una distribución $\mathsf{N}(\mu, \sigma)$ contra +los cuantiles de una distribución $\mathsf{N}(0,1)$, estos cuantiles aparecen en una línea recta: + +```{r, fig.width=4,fig.height=3} +comparacion_tbl <- tibble(f = seq(0.01, 0.99, 0.01)) |> + mutate(cuantiles_normal = qnorm(f, 5, 3), + cuantiles_norm_estandar = qnorm(f, 0, 1)) +ggplot(comparacion_tbl, aes(cuantiles_norm_estandar, cuantiles_normal)) + + geom_point() +``` + +Ahora supongamos que tenemos una muestra $X_1, \ldots, X_n$. ¿Cómo podemos +checar si estos datos tienen una distribución aproximadamente normal? + +- Si la muestra tiene una distribución aproximadamente $\mathsf{N}(\mu, \sigma)$, entonces +sus cuantiles muestrales y los cuantiles respectivos de la normal estándar están aproximadamente +en una línea recta. + +Primero veamos un ejemplo donde los datos son generados según una normal. + +```{r, fig.width=7, fig.height=3} +set.seed(21) +muestra <- tibble(x_1 = rnorm(60, 10, 3), x_2 = rgamma(60, 2, 5)) +graf_1 <- ggplot(muestra, aes(sample = x_1)) + + geom_qq(distribution = stats::qnorm) + + geom_qq_line(colour = "red") +graf_2 <- ggplot(muestra, aes(sample = x_2)) + + geom_qq(distribution = stats::qnorm) + + geom_qq_line(colour = "red") +graf_1 + graf_2 +``` +¿Cuáles son los datos aproximadamente normales? ¿Cómo interpretas las desviaciones de +la segunda gráfica en términos de la forma de la distribución normal? + +## Prueba de hipótesis de normalidad {-} + +Para interpretar las gráficas de cuantiles normales se requiere práctica, +pues claramente los datos, aún cuando provengan de una distribución normal, no +van a caer justo sobre una línea recta y observaremos variabilidad. Esto no descarta +necesariamente que los datos sean aproximadamente normales. Con la práctica, generalmente +esta gráfica nos da una buena indicación si el supuesto de normalidad es apropiado. + +Sin embargo, podemos hacer una prueba de hipótesis formal de normalidad si +quisiéramos. La hipótesis nula es la siguiente: + +- Los datos provienen de una distribución normal, y las desviaciones que observamos +de una línea recta se deben a variación muestral. +- Podemos generar datos nulos tomando la media y desviación estándar muestrales, y +generando muestras normales $\mathsf{N}(\bar{x}, s)$. +- Usamos el *lineup*, produciendo datos bajo la hipótesis nula y viendo si +podemos distinguir los datos. Por ejemplo: + +```{r, fig.height=7.2, fig.width=7} +library(nullabor) +lineup_normal <- lineup(null_dist("x_2", dist = "normal"), muestra) +ggplot(lineup_normal, aes(sample = x_2)) + + geom_qq(distribution = stats::qnorm) + + geom_qq_line(colour = "red") + + facet_wrap(~ .sample) +``` + +En esta gráfica claramente rechazaríamos la hipótesis de normalidad. Sin embargo, para +la primera muestra, obtenemos: + +```{r, fig.height=7.2, fig.width=7} +lineup_normal <- lineup(null_dist("x_1", dist = "normal"), muestra) +ggplot(lineup_normal, aes(sample = x_1)) + + geom_qq(distribution = stats::qnorm) + + geom_qq_line(colour = "red") + + facet_wrap(~ .sample) +``` +Los datos verdaderos están en + +```{r} +attr(lineup_normal, "pos") +``` + +### Ejemplo {-} + +Consideremos el problema de estimar el total poblacional de los precios de las casas +que se vendieron. El estimador que usamos fue la suma muestral expandida por un +factor. Vamos a checar qué tan cerca de la normalidad está la distribución de +meustreo de esta estadística ($n=250$): + +```{r} +replicaciones_2 +``` + + +```{r, fig.width = 4, fig.height=3} +ggplot(replicaciones_2, aes(sample = est_total_millones)) + + geom_qq(alpha = 0.3) + + geom_qq_line(colour = "red") +``` +Y vemos que en efecto el TCL aplica en este ejemplo, y la aproximación es buena. +Aunque la población original es sesgada, la descripción de la distribución de +muestreo es sorprendemente compacta: + +- La distribución de muestreo de nuestro estimador del total $\hat{t}$ es +aproximadamente normal con media $\bar{x}$ y desviación estándar $s$, donde: + +```{r} +mu <- mean(replicaciones_2$est_total_millones) +s <- sd(replicaciones_2$est_total_millones) +c(mu = mu, s = s) |> round(2) +``` +Estas cantidades están en millones de dólares. + +### Ejemplo {-} + +Supongamos que queremos calcular la probabilidad que la suma de 30 variables +uniformes en $[0,1]$ independientes sea mayor que 18. Podríamos aproximar esta +cantidad usando simulación. Otra manera de aproximar esta cantidad es con +el TCL, de la siguiente forma: + +Si $S=X_1 + X_2 + X_{30}$, entonces la media de $S$ es 15 (¿cómo se calcula?) +y su desviación estándar es $\sqrt{\frac{30}{12}}$. La suma es entonces +aproximadamente $\mathsf{N}\left(15, \sqrt{\frac{30}{12}}\right)$. Entonces + +$$P(S > 18) = P \left (\frac{S - 15}{\sqrt{\frac{30}{12}}} > \frac{18 - 15}{\sqrt{\frac{30}{12}}}\right) \approx P(Z > 1.897)$$ + +donde $Z$ es normal estándar. Esta última cantidad la calculamos +usando la función de distribución de la normal estándar, y nuestra aproximación es + +```{r} +1 - pnorm(1.897) +``` + +Podemos checar nuestro cálculo usando simulación: + +```{r} +tibble(n_sim = 1:100000) |> + mutate(suma = map_dbl(n_sim, ~ sum(runif(30)))) |> + summarise(prob_may_18 = mean(suma > 18), .groups = "drop") +``` +Y vemos que la aproximación normal es buena para fines prácticos. + +```{block, type='ejercicio'} +Usando simulaciones haz un histograma que aproxime la distribución de muestreo de $S$. +Haz una gráfica de cuantiles normales para checar la normalidad de esta distribución. +``` + +## Ejemplo {-} +Cuando el sesgo de la distribución poblacional es grande, puede ser necesario +que $n$ sea muy grande para que la aproximación normal sea aceptable para el +promedio o la suma. Por ejemplo, si tomamos una gamma con parámetro de forma chico, +$n = 30$ no es suficientemente bueno, especialmente si quisiéramos +aproximar probabilidades en las colas de la distribución: + + +```{r} +sims_gamma <- map_df(1:2000, ~ tibble(suma = mean(rgamma(30, 0.1, 1))), + .id = "n_sim") +ggplot(sims_gamma, aes(x = suma)) + geom_histogram() +``` + +## Más del Teorema central del límite {-} + +- El teorema central del límite aplica a situaciones más generales que +las del enunciado del teorema básico. Por ejemplo, aplica a poblaciones +finitas (como vimos en el ejemplo de las casas) bajo muestreo sin +reemplazo, y aplica también a otras estadísticas como los cuantiles muestrales. + +- Es importante notar que la calidad de la aproximación del TCL depende de características +de la población y también del tamaño de muestra $n$. Para ver si el TCL aplica, podemos hacer ejercicios de simulación bajo diferentes supuestos acerca de la población. +También veremos más adelante, con remuestreo, maneras de checar si es factible el TCL dependiendo del análisis de una muestra dada que tengamos. + +- El TCL era particularmente importante en la práctica antes de que pudiéramos +hacer simulación por computadora. Era la única manera de aproximar y entender la distribución muestral fuera de cálculos analíticos (como los que hicimos para +el máximo de un conjunto de uniformes, por ejemplo). + +- Hoy en día, veremos que podemos hacer simulación para obtener respuestas más +exactas, particularmente en la construcción de intervalos de confianza, por ejemplo. Dependemos menos de **resultados asintóticos**, como el TCL. + +- Cuando aproximamos una distribución discreta mediante la distribución normal, +conviene hacer *correcciones de continuidad*, como se explica en [@Chihara], 4.3.2. + + + + + + + diff --git a/89-transformaciones.Rmd b/89-transformaciones.Rmd index 6570da6..0933a0d 100644 --- a/89-transformaciones.Rmd +++ b/89-transformaciones.Rmd @@ -12,30 +12,29 @@ library(patchwork) En ocasiones es conveniente transformar los datos para el análisis, el objetivo de los ajustes es simplificar la interpretación y el análisis al -eliminar fuentes de variación conocidas, o para simplificar los patrones. +eliminar fuentes de variación conocidas, también es común realizan +transformaciones para simplificar los patrones. Algunos ejemplos donde eliminamos efectos conocidos: 1. Cuando analizamos el precio de venta de las casas podemos eliminar la variación debida al tamaño de las casas al pasar de precio de venta a precio de venta por metro cuadrado. -De manera similar considerar propina como porcentaje de la cuenta. +De manera similar al analizar las propinas puede convenir considerar la propina como porcentaje de la cuenta. 2. En series de tiempo cuando los datos están relacionados con el tamaño de la población podemos ajustar a mediciones per capita (en series de tiempo PIB). También es común ajustar por inflación, o poner cantidades monetarias en valor presente. -```{r} +```{r, fig.height=8, fig.width=6} mex_dat <- global_economy |> filter(Code == "MEX") pib <- ggplot(mex_dat, aes(x = Year, y = GDP / 1e6)) + geom_line() -``` -```{r, out.height=4, out.width=7.5} pib_pc <- ggplot(mex_dat, aes(x = Year, y = GDP / Population)) + geom_line() -pib + pib_pc +pib / pib_pc ``` @@ -44,23 +43,32 @@ el patrón en los datos y la interpretación. Veamos un ejemplo donde es apropiado la transformación logaritmo. -```{r} +Usamos los datos Animals con información de peso corporal promedio y peso cerebral +promedio para 28 especies. Buscamos entender la relación entre estas dos +variables, e inspeccionar que especies se desvían (residuales) del esperado. +Comenzamos con un diagrama de dispersión usando las unidades originales + +```{r, warning=FALSE, fig.height=5} animals_tbl <- as_tibble(Animals, rownames = "animal") p1 <- ggplot(animals_tbl, aes(x = body, y = brain, label = animal)) + - geom_point() + geom_point() + + labs(subtitle = "Unidades originales") p2 <- ggplot(animals_tbl, aes(x = body, y = brain, label = animal)) + geom_point() + xlim(0, 500) + ylim(0, 1500) + - geom_text_repel() + geom_text_repel() + + labs(subtitle = "Unidades originales, eliminando 'grandes'") (p1 + p2) ``` -Cuando hacemos la transformación logaritmo obtenemos una gráfica más fácil de -leer y los datos se modelarán con más facilidad. +Incluso cuando nos limitamos a especies de menos de 500 kg de masa corporal, la +relación no es fácil de descrubir.En la suguiente gráfica hacemos la transformación +logaritmo y obtenemos una gráfica más fácil de +leer, además los datos se modelarán con más facilidad. -```{r} +```{r, warning=FALSE} p3 <- ggplot(animals_tbl, aes(x = log(body), y = log(brain), label = animal)) + geom_smooth(method = "lm", se = FALSE, color = "red") + geom_point() + @@ -70,7 +78,7 @@ p3 <- ggplot(animals_tbl, aes(x = log(body), y = log(brain), label = animal)) + p3 ``` -La transformación logaritmo tiene además ventajas en interpretación, para diferencias +La transformación logaritmo tiene también ventajas en interpretación, para diferencias chicas en escala log, las diferencias corresponden a diferencias porcentuales en la escala original, por ejempo consideremos la diferencia entre el peso en escala log de humano y borrego: 4.13 - 4.02 = 0.11. @@ -79,9 +87,11 @@ escala original: 62/55.5 - 1 = 0.12 ```{r} -animals_tbl |> +animals_tbl <- animals_tbl |> mutate(log_body = log(body), - log_brain = log(brain)) |> + log_brain = log(brain)) + +animals_tbl |> filter(animal == "Human" | animal == "Sheep") |> arrange(body) |> gt::gt() |> @@ -89,17 +99,32 @@ animals_tbl |> ``` + +```{r, eval=FALSE, include=FALSE, echo=FALSE} +fit <- lm(log_brain ~ log_body, data = animals_tbl) +animals_tbl |> + mutate(log_body = log(body), + log_brain = log(brain), + log_brain_pred = fitted(fit)) |> + arrange(body) +``` + + + Y podemos usarlo también para interpretar la recta de referencia $y = 2.55 + 0.5 x$ , para cambios chicos: *Un incremento de 10% en masa total corresponde en un incremento de 5% en masa cerebral.* +El coeficiente de la regresión log-log, en nuestro ejemplo 0.5, es la [elasticidad](https://en.wikipedia.org/wiki/Elasticity_(economics)) y es un concepto +común en economía. + **Justificación** Para entender la interpretación como cambio porcentual recordemos primero que la representación con series de Taylor de la función exponencial es: -$$e^x = \sum_{n=1}^\infty \frac{x^n}{n!}$$ +$$e^x = \sum_{n=0}^\infty \frac{x^n}{n!}$$ Más aún podemos tener una aproximación usando [polinomios de Taylor](https://en.wikipedia.org/wiki/Taylor%27s_theorem), en el caso de la exponencial el $k$-ésimo polinomio de Taylor está dado por: @@ -111,8 +136,8 @@ razonable y tenemos: $$Ae^{\delta} \approx A(1+\delta)$$ -```{r} -dat <- tibble(delta = seq(0, 0.25, 0.01), exp_delta = exp(delta), uno_mas_delta = 1 + delta) +```{r, fig.height=4, fig.width=4.5} +dat <- tibble(delta = seq(0, 1, 0.01), exp_delta = exp(delta), uno_mas_delta = 1 + delta) ggplot(dat, aes(x = uno_mas_delta, y = exp_delta)) + geom_line() + diff --git a/renv.lock b/renv.lock index cc81170..cf022e5 100644 --- a/renv.lock +++ b/renv.lock @@ -9,6 +9,13 @@ ] }, "Packages": { + "BH": { + "Package": "BH", + "Version": "1.84.0-0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "a8235afbcd6316e6e91433ea47661013" + }, "DBI": { "Package": "DBI", "Version": "1.2.3", @@ -73,6 +80,19 @@ ], "Hash": "1920b2f11133b12350024297d8a4ff4a" }, + "MatrixModels": { + "Package": "MatrixModels", + "Version": "0.5-3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Matrix", + "R", + "methods", + "stats" + ], + "Hash": "0776bf7526869e0286b0463cb72fb211" + }, "R6": { "Package": "R6", "Version": "2.5.1", @@ -104,6 +124,32 @@ ], "Hash": "f27411eb6d9c3dada5edd444b8416675" }, + "SparseM": { + "Package": "SparseM", + "Version": "1.84-2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "graphics", + "methods", + "stats", + "utils" + ], + "Hash": "e78499cbcbbca98200254bd171379165" + }, + "TH.data": { + "Package": "TH.data", + "Version": "1.1-2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "MASS", + "R", + "survival" + ], + "Hash": "5b250ad4c5863ee4a68e280fcb0a3600" + }, "V8": { "Package": "V8", "Version": "4.4.2", @@ -129,6 +175,18 @@ ], "Hash": "4f57884290cc75ab22f4af9e9d4ca862" }, + "anytime": { + "Package": "anytime", + "Version": "0.3.9", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "BH", + "R", + "Rcpp" + ], + "Hash": "74a64813f17b492da9c6afda6b128e3d" + }, "askpass": { "Package": "askpass", "Version": "1.2.0", @@ -229,6 +287,18 @@ ], "Hash": "896a79478a50c78fb035a37148638f4e" }, + "boot": { + "Package": "boot", + "Version": "1.3-30", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "graphics", + "stats" + ], + "Hash": "96abeed416a286d4a0f52e550b612343" + }, "broom": { "Package": "broom", "Version": "1.0.6", @@ -406,6 +476,16 @@ ], "Hash": "1b4ced11b3b6c23f0e90a6becee4d983" }, + "codetools": { + "Package": "codetools", + "Version": "0.2-20", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "61e097f35917d342622f21cdc79c256e" + }, "colorspace": { "Package": "colorspace", "Version": "2.1-1", @@ -427,6 +507,18 @@ "Repository": "CRAN", "Hash": "5d8225445acb167abf7797de48b2ee3c" }, + "confintr": { + "Package": "confintr", + "Version": "1.0.2", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "boot", + "stats" + ], + "Hash": "f3e7567eabebba308f9c825a2bf249cf" + }, "conflicted": { "Package": "conflicted", "Version": "1.2.0", @@ -607,6 +699,17 @@ ], "Hash": "4ef372b716824753719a8a38b258442d" }, + "ellipsis": { + "Package": "ellipsis", + "Version": "0.3.2", + "Source": "Repository", + "Repository": "RSPM", + "Requirements": [ + "R", + "rlang" + ], + "Hash": "bb0eec2fe32e88d9e2836c2f73ea2077" + }, "evaluate": { "Package": "evaluate", "Version": "0.24.0", @@ -822,6 +925,63 @@ ], "Hash": "44c6a2f8202d5b7e878ea274b1092426" }, + "ggpmisc": { + "Package": "ggpmisc", + "Version": "0.6.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "MASS", + "R", + "confintr", + "dplyr", + "generics", + "ggplot2", + "ggpp", + "grid", + "lmodel2", + "lubridate", + "multcomp", + "multcompView", + "plyr", + "polynom", + "quantreg", + "rlang", + "scales", + "splus2R", + "stats", + "tibble" + ], + "Hash": "e0e669a53c9e39482705df3de75ec054" + }, + "ggpp": { + "Package": "ggpp", + "Version": "0.5.8-1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "MASS", + "R", + "dplyr", + "ggplot2", + "glue", + "grDevices", + "grid", + "gridExtra", + "lubridate", + "magrittr", + "polynom", + "rlang", + "scales", + "stats", + "stringr", + "tibble", + "vctrs", + "xts", + "zoo" + ], + "Hash": "85f25ea0bb52935bd33d30b2c6c8f0a6" + }, "ggrepel": { "Package": "ggrepel", "Version": "0.9.5", @@ -921,6 +1081,20 @@ ], "Hash": "d6db1667059d027da730decdc214b959" }, + "gridExtra": { + "Package": "gridExtra", + "Version": "2.3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "grDevices", + "graphics", + "grid", + "gtable", + "utils" + ], + "Hash": "7d7f283939f563670a697165b2cf5560" + }, "gt": { "Package": "gt", "Version": "0.11.0", @@ -1214,6 +1388,16 @@ ], "Hash": "b8552d117e1b808b09a832f589b79035" }, + "lmodel2": { + "Package": "lmodel2", + "Version": "1.7-3", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R" + ], + "Hash": "d6402cf7948749ef56bc8d9c8099bfbe" + }, "lpSolve": { "Package": "lpSolve", "Version": "5.6.20", @@ -1356,6 +1540,32 @@ "Repository": "CRAN", "Hash": "bb94fd8ee5f7f127eae2bddbff26864d" }, + "multcomp": { + "Package": "multcomp", + "Version": "1.4-26", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "TH.data", + "codetools", + "graphics", + "mvtnorm", + "sandwich", + "stats", + "survival" + ], + "Hash": "ec6f951a557132215fab91912acdd9ef" + }, + "multcompView": { + "Package": "multcompView", + "Version": "0.1-10", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "grid" + ], + "Hash": "bd174e328b679e9a1ce589422002c72b" + }, "munsell": { "Package": "munsell", "Version": "0.5.1", @@ -1367,6 +1577,17 @@ ], "Hash": "4fd8900853b746af55b81fda99da7695" }, + "mvtnorm": { + "Package": "mvtnorm", + "Version": "1.2-6", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "stats" + ], + "Hash": "b5b11e71e41d9d0eef7d8dc26903ed55" + }, "nlme": { "Package": "nlme", "Version": "3.1-165", @@ -1476,6 +1697,28 @@ ], "Hash": "01f28d4278f15c76cddbea05899c5d6f" }, + "plyr": { + "Package": "plyr", + "Version": "1.8.9", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "Rcpp" + ], + "Hash": "6b8177fd19982f0020743fadbfdbd933" + }, + "polynom": { + "Package": "polynom", + "Version": "1.4-1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "graphics", + "stats" + ], + "Hash": "ceb5c2a59ba33d42d051285a3e8a5118" + }, "posterior": { "Package": "posterior", "Version": "1.6.0", @@ -1585,6 +1828,24 @@ ], "Hash": "1cba04a4e9414bdefc9dcaa99649a8dc" }, + "quantreg": { + "Package": "quantreg", + "Version": "5.98", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "MASS", + "Matrix", + "MatrixModels", + "R", + "SparseM", + "graphics", + "methods", + "stats", + "survival" + ], + "Hash": "017561f17632c065388b7062da030952" + }, "ragg": { "Package": "ragg", "Version": "1.3.2", @@ -1815,6 +2076,19 @@ ], "Hash": "3c8013cdd7f1d20de5762e3f97e5e274" }, + "sandwich": { + "Package": "sandwich", + "Version": "3.1-0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "stats", + "utils", + "zoo" + ], + "Hash": "1cf6ae532f0179350862fefeb0987c9b" + }, "sass": { "Package": "sass", "Version": "0.4.9", @@ -1885,6 +2159,17 @@ ], "Hash": "ad57b543f7c3fca05213ba78ff63df9b" }, + "splus2R": { + "Package": "splus2R", + "Version": "1.3-5", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "methods" + ], + "Hash": "51f4ca9989e0de9782037e426298bee6" + }, "stringi": { "Package": "stringi", "Version": "1.8.4", @@ -1915,6 +2200,22 @@ ], "Hash": "960e2ae9e09656611e0b8214ad543207" }, + "survival": { + "Package": "survival", + "Version": "3.7-0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "Matrix", + "R", + "graphics", + "methods", + "splines", + "stats", + "utils" + ], + "Hash": "5aaa9cbaf4aba20f8e06fdea1850a398" + }, "svglite": { "Package": "svglite", "Version": "2.1.3", @@ -2118,6 +2419,41 @@ ], "Hash": "584eda9966ed571b3adc75e338116147" }, + "tsibble": { + "Package": "tsibble", + "Version": "1.1.5", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "anytime", + "dplyr", + "ellipsis", + "generics", + "lifecycle", + "lubridate", + "methods", + "rlang", + "tibble", + "tidyselect", + "vctrs" + ], + "Hash": "a75e397766b45996310908b5b32557ba" + }, + "tsibbledata": { + "Package": "tsibbledata", + "Version": "0.4.1", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "rappdirs", + "tsibble", + "utils", + "vctrs" + ], + "Hash": "6b53c37a18049320fb08478aef6310c2" + }, "tweenr": { "Package": "tweenr", "Version": "2.0.3", @@ -2272,12 +2608,39 @@ ], "Hash": "1d0336142f4cd25d8d23cd3ba7a8fb61" }, + "xts": { + "Package": "xts", + "Version": "0.14.0", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "methods", + "zoo" + ], + "Hash": "288bc0fd02bb4e1489f37831d8803b2b" + }, "yaml": { "Package": "yaml", "Version": "2.3.10", "Source": "Repository", "Repository": "CRAN", "Hash": "51dab85c6c98e50a18d7551e9d49f76c" + }, + "zoo": { + "Package": "zoo", + "Version": "1.8-12", + "Source": "Repository", + "Repository": "CRAN", + "Requirements": [ + "R", + "grDevices", + "graphics", + "lattice", + "stats", + "utils" + ], + "Hash": "5c715954112b45499fb1dadc6ee6ee3e" } } }