Loading 'brms' package (version 2.22.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
Attaching package: 'brms'
The following object is masked from 'package:stats':
ar
library(CalvinBayes)
Loading required package: dplyr
Attaching package: 'dplyr'
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
Loading required package: ggplot2
Loading required package: ggformula
Loading required package: scales
Loading required package: ggridges
New to ggformula? Try the tutorials:
learnr::run_tutorial("introduction", package = "ggformula")
learnr::run_tutorial("refining", package = "ggformula")
Loading required package: bayesplot
This is bayesplot version 1.11.1
- Online documentation and vignettes at mc-stan.org/bayesplot
- bayesplot theme set to bayesplot::theme_default()
* Does _not_ affect other ggplot2 plots
* See ?bayesplot_theme_set for details on theme setting
Attaching package: 'bayesplot'
The following object is masked from 'package:brms':
rhat
Attaching package: 'CalvinBayes'
The following object is masked from 'package:bayesplot':
rhat
The following objects are masked from 'package:brms':
rhat, rstudent_t
The following object is masked from 'package:datasets':
HairEyeColor
library(coda)library(mosaic)
Registered S3 method overwritten by 'mosaic':
method from
fortify.SpatialPolygonsDataFrame ggplot2
The 'mosaic' package masks several functions from core packages in order to add
additional features. The original behavior of these functions should not be affected by this.
Attaching package: 'mosaic'
The following object is masked from 'package:Matrix':
mean
The following object is masked from 'package:scales':
rescale
The following object is masked from 'package:ggplot2':
stat
The following objects are masked from 'package:dplyr':
count, do, tally
The following object is masked from 'package:brms':
mm
The following objects are masked from 'package:stats':
binom.test, cor, cor.test, cov, fivenum, IQR, median, prop.test,
quantile, sd, t.test, var
The following objects are masked from 'package:base':
max, mean, min, prod, range, sample, sum
options(width =100)set.seed(100)
9.1 Modelo Logístico
Supongamos que queremos predecir el género de una persona a partir de su altura y peso. ¿Qué tipo de modelo deberíamos usar? Dada una combinación particular de altura y peso, una persona puede ser hombre o mujer. Por lo tanto, para una combinación dada de altura/peso, la distribución del género es una variable aleatoria Bernoulli. Nuestro objetivo es convertir la combinación de altura/peso en el parámetro \(\theta\), que especifica la proporción de personas con esa altura y peso que son hombres (o mujeres).
Entonces, nuestro modelo se ve algo así:
\[\begin{align*}
Y &\sim \text{Bernoulli}(\theta) \\
\theta &= \mathrm{f}(\texttt{altura}, \texttt{peso})
\end{align*}\]
Pero, ¿qué funciones deberíamos usar para \(f\)?
9.1.1 El enfoque habitual: Regresión logística
Nos gustaría usar una función lineal, pero convertir su rango de \((-\infty, \infty)\) a \((0, 1)\). Alternativamente, podemos convertir el rango \((0,1)\) a \((-\infty, \infty)\) y parametrizar la distribución Bernoulli de manera diferente.
La transformación más común utiliza la transformación de logaritmo de odds:
La función logit se llama función de enlace y la función logística es la función de enlace inversa.
9.1.2 Otros enfoques
Podríamos hacer algo similar con cualquier par de funciones que conviertan hacia adelante y hacia atrás entre \((0,1)\) y \((-\infty, \infty)\). Para cualquier variable aleatoria, la función de distribución acumulativa (CDF) tiene dominio \((-\infty, \infty)\) y rango \((0,1)\), por lo que
Cualquier par de funciones CDF/inversa de CDF puede ser usado en lugar de la transformación logística.
Usar pnorm() y qnorm() se llama regresión probit.
gf_function(ilogit, xlim =c(-6, 6), color =~"logística") %>%gf_function(pnorm, color =~"probit (estándar)",xlim =c(-6, 6)) %>%gf_function(pnorm, args =list(sd =1.8), color =~"probit (media = 0, sd = 1.8)",xlim =c(-6, 6)) %>%gf_theme(legend.position ="top") %>%gf_labs(color ="")
9.1.3 Preparación de los datos
Utilizaremos un subconjunto de los datos de NHANES para este ejemplo. Dado que los datos están en libras y pulgadas, convertiremos los datos de NHANES a estas unidades.
Como en otros modelos, necesitamos convertir nuestra variable dicotómica en 0’s y 1’s. Sin duda, queremos excluir a los niños del modelo, ya que los patrones de altura y peso son diferentes para niños y adultos. De hecho, seleccionaremos únicamente a los individuos de 22 años.
Aprovecharemos también para eliminar las variables que no necesitamos y eliminar cualquier fila que tenga valores faltantes en esas tres columnas.
library(NHANES)library(brms)nhanes <- NHANES %>%mutate(weight = Weight *2.2,height = Height /2.54,male =as.numeric(Gender =="male") ) %>%filter(Age ==22) %>%select(male, height, weight) %>%filter(complete.cases(.)) # eliminar filas con valores faltantes en alguna de las 3 variables
9.1.4 Especificación de la familia y función de enlace en brm()
Comparado con nuestro modelo habitual de regresión lineal, necesitamos hacer dos ajustes:
Utilizar la familia de distribuciones Bernoulli para el ruido.
Utilizar la función de enlace logit (inversa logística) para traducir hacia adelante y hacia atrás entre la parte lineal del modelo y la distribución.
Así que para la regresión logística y la regresión probit, utilizamos:
logistic_brm <-brm(male ~ height + weight, family =bernoulli(link = logit), data = nhanes)
probit_brm <-brm(male ~ height + weight, family =bernoulli(link = probit), data = nhanes)
El resto del modelo se comporta como antes. En particular, se utiliza una distribución t para la intersección y distribuciones uniformes impropias para los demás coeficientes de regresión. Este modelo no tiene un parámetro \(\sigma\). (La varianza de una distribución Bernoulli está determinada por el parámetro de probabilidad.):
prior_summary(logistic_brm)
9.1.5 Interpretación de Modelos de Regresión Logística
Antes de analizar nuestro modelo con dos predictores, veamos un modelo con solo un predictor.
male_by_weight_brm <-brm(male ~ weight, family =bernoulli(), data = nhanes)
Claramente hay una tendencia ascendente: las personas más pesadas tienen más probabilidades de ser hombres que las personas más ligeras. Esto se observa en la distribución posterior para \(\beta_{\mathrm{weight}}\).
hdi(posterior(male_by_weight_brm), pars ="b_weight")
El conjunto de datos de pulpa del paquete faraway consiste en información sobre el brillo de la pulpa de papel según el operador de turno a-d. Las variables incluidas son brillo y operador. En este contexto, un operador de turno se refiere a un trabajador que opera una máquina o supervisa un proceso durante un turno específico en una planta de producción de pulpa de papel. Los turnos están generalmente divididos en diferentes períodos del día, y los operadores de turno a, b, c y d se encargan de las operaciones durante estos diferentes períodos. El desempeño de estos operadores puede afectar variables como el brillo de la pulpa de papel producida.
Ilustramos con un experimento para comparar cuatro procesos, A, B, C y D, para la producción de penicilina. Estos son los tratamientos. La materia prima, licor de maíz, es bastante variable y solo puede ser preparada en mezclas suficientes para cuatro corridas. Por lo tanto, un diseño de bloques completos aleatorizados es sugerido por la naturaleza de las unidades experimentales.
data("penicillin")penicillin$Blend <-gl(5,4)
Gráfico comparando la producción media de penicilina con los otros dos factores: