Show code
script src = "https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js"Show code
$(document).ready(function() {
$('body').prepend('<div class=\"zoomDiv\"><img src=\"\" class=\"zoomImg\"></div>');
// onClick function for all plots (img's)
$('img:not(.zoomImg)').click(function() {
$('.zoomImg').attr('src', $(this).attr('src')).css({width: '100%'});
$('.zoomDiv').css({opacity: '1', width: 'auto', border: '1px solid white', borderRadius: '5px', position: 'fixed', top: '50%', left: '50%', marginRight: '-50%', transform: 'translate(-50%, -50%)', boxShadow: '0px 0px 50px #888888', zIndex: '50', overflow: 'auto', maxHeight: '100%'});
});
// onClick function for zoomImg
$('img.zoomImg').click(function() {
$('.zoomDiv').css({opacity: '0', width: '0%'});
});
});Show code
knitr::opts_chunk$set(warning = TRUE, message = TRUE)
rm(list=ls());gc()
used (Mb) gc trigger (Mb) max used (Mb)
Ncells 646813 34.6 1332658 71.2 904299 48.3
Vcells 1170604 9.0 8388608 64.0 2381499 18.2Show code
knitr::opts_chunk$set(
echo = TRUE,
message = TRUE,
warning = TRUE
)
local({r <- getOption("repos")
r["CRAN"] <- "http://cran.r-project.org"
options(repos=r)
})
if(!require(pacman)){install.packages("pacman")}
Loading required package: pacmanWarning: package 'pacman' was built under R version 4.1.3Show code
pacman::p_load(devtools, here, showtext, ggpattern, RefManageR, pagedown, magick, bibtex, DiagrammeR, xaringan, xaringanExtra, xaringanthemer, fontawesome, widgetframe, datapasta, tidyverse, psych, cowplot, coxphw, future, timereg, flexsurv, pdftools, mstate, showtext, compareGroups, chilemapas, choroplethrAdmin1, choroplethr, choroplethrMaps, ggiraph, sf,distill, qrcode, pdftools, dagitty, ggdag, geepack, survey, ipw, reshape, widgetframe, geeM, MuMIn, confoundr, install=F)
copiar_nombres <- function(x,row.names=FALSE,col.names=TRUE,dec=",",...) {
if(class(ungroup(x))[1]=="tbl_df"){
if(options()$OutDec=="."){
options(OutDec = dec)
write.table(format(data.frame(x)),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
options(OutDec = ".")
return(x)
} else {
options(OutDec = ",")
write.table(format(data.frame(x)),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
options(OutDec = ",")
return(x)
}
} else {
if(options()$OutDec=="."){
options(OutDec = dec)
write.table(format(x),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
options(OutDec = ".")
return(x)
} else {
options(OutDec = ",")
write.table(format(x),"clipboard",sep="\t",row.names=FALSE,col.names=col.names,...)
options(OutDec = ",")
return(x)
}
}
}
options(knitr.kable.NA = '')
Introducción
Causalidad
Muchas preguntas en salud pública se remiten a preguntas causales. Desde determinantes sociales (por ej., origen de las inequidades, vulnerabilidad social, políticas públicas), pasando por factores de riesgo, procesos biológicos/mecanismos etiológicos, etc. Estas preguntas no se limitan a describir sino que algunas apuntan a un análisis causal.1 Al margen de la tarea de la vigilancia involucrada en el ejercicio descriptivo, en algunos casos las preguntas requieren imaginar un escenario hipotético a futuro, o ¿qué pasaría si hacemos/aplicamos \(A\) política pública? (dirigida a cuantificar los efectos de causas), vs. inferencia al pasado, o ¿qué causa \(Y\)? (dirigida a identificar las causas de los efectos)”.2
Además de que tales preguntas refieren a escenarios que no pudimos ni podríamos ver directamente en la práctica, también establecen vínculos entre uno o más elementos. Pero no una relación de cualquier tipo, sino uno de carácter causal. Hay asociación cuando el tener información sobre una variable \(A\) me permite predecir mejor en promedio el resultado de otra, incluso si esa variable no causa \(Y\). El mecanismo no lo sabemos. Uno de los desafíos es poder determinar si por sí sola una variable \(A\) causa \(Y\), ya que muchas veces hay terceras variables (\(L\)) que podrían explicar las asociaciones entre \(A\) e \(Y\).3,4
Estos métodos de identificación causal se enmarcan en una aproximación de causalidad conocida como “potential outcomes framework”.5 Este marco de inferencia contrafactual iniciado por Neyman y continuado por Rubin y Robins, busca comparar teóricamente qué ocurriría si todos los pacientes hubieran sido tratados frente al caso que ningún paciente hubiese sido tratado6. Un diseño de investigación casi-perfecto para este marco de investigación son los ensayos clínicos controlados, ya que al ser aleatorizados, en la mayoría de los casos permiten estimar el impacto de un determinado tratamiento o intervención por sí misma no sólo en los tratados, si no también el efecto que hubiera tenido en los no tratados de haberlo sido.7,8
No obstante, en muchos casos no es posible, ya sea por razones prácticas (factibilidad económica o social, por ejemplo, para lograr potencia estadística, o en el caso del efecto de un shock económico o una catástrofe natural, el tratamiento no es del todo manipulable) o ética (por ejemplo, no puedo inducir a que sujetos consuman un producto del que se sospecha un perjuicio para la salud). De igual manera, en algunos casos también los estudios aleatorizados tendrán problemas para estimar efectos en caso de pérdida de seguimiento, variables perdidas, falta de cooperación, entre otros. En tal escenario se hace necesario recurrir a estudios observacionales.6,9 El control de la confusión es uno de los pasos más desafiantes en la conducción de estudios observacionales y es aún más importante el cumplimiento de ciertas condiciones para dar validez a nuestras inferencias causales.10
El planteamiento de Pearl presentado a continuación permite interpretar que los datos de por sí no son suficientes para señalar que una relación es causal. Como diría Hernán6, la interpretación depende de consideraciones más allá de los datos, como por ejemplo, el conocimiento adquirido:
“Estas consideraciones implican que el eslogan ‘correlación no implica causalidad’ puede traducirse en un principio útil: uno no puede generar declaraciones causales sólo por asociaciones, incluso con datos poblacionales. Cada conclusión causal debe descansar detrás de supuestos causales no medibles en estudios observacionales” [These considerations imply that the slogan “correlation does not imply causation” can be translated into a useful principle: one cannot substantiate causal claims from associations alone, even at the population level—behind every causal conclusion there must lie some causal assumption that is not testable in observational studies]11 (p.4)
Supuestos
Uno de los principales supuestos son acuñados por Rubin en 1980 y al que denomina “la hipótesis de estabilidad del valor del tratamiento unitario” o stable-unit-treatment-value (SUTVA por sus siglas en inglés).12 La validez de nuestras inferencias causales requieren el cumplimiento de dichas condiciones13. Lo anterior nos lleva a pensar el mecanismo de asignación detrás de los datos obtenidos (assignment mechanism) para definir efectos causales a identificar. Si la identificación es fallida, no se puede aprender una relación causal de los datos incluso teniendo una muestra infinita14. A continuación se describen:3,13,15
Consistencia: Constituye uno de los criterios más cercanos a SUTVA Si un individuo tiene una exposición \(A=a\), entonces su resultado \(Y = Y(a)\). Dicho de otra manera, que la misma exposición tiene el mismo efecto en un individuo sin importar cómo llegó a ser expuesto. Implica hacerse las siguientes preguntas: ¿Todos recibieron el tratamiento en medida similar?, ¿Lo mismo ocurre con quienes no lo recibieron?, ¿Está bien clasificada la exposición?, ¿Es ambigua la definición de la exposición?. Nos permite homologar los resultados de quienes recibieron una exposición al resultado potencial que se habría observado con ese valor de exposición. Por otra parte, se requiere que no existiera interferencia (que el resultado potencial no estuviera afectado por la exposición a tratamiento por otros sujetos) o que el efecto de la exposición en un individuo sea independiente de la exposición en otros. De igual manera, se debe asumir que no hay variaciones ocultas del tratamiento o versiones ocultas que puedan tener otro efecto9,12,16
Positividad: La identificación de efectos causales requieren una variabilidad suficiente en la exposición o en la asignación de tratamiento en los estratos de los confusores. Lo anterior implica hacerse las siguientes preguntas: ¿La distribución de las covariables \(L\) en tratados es similar que los no-tratados?, ¿cada observación tiene el potencial de estar (o no) expuesto a cada nivel de intervención?. La probabilidad de asignación está comprendida entre 0 y 1, por lo que se puede tener tratados y no tratados. Es necesario que haya expuestos y no expuestos en todos los niveles de confusión. En otras palabras, no puede haber un confusor o covariable que determine perfectamente el tratamiento o no-tratamiento fuera de \(A\).3,16,17
Intercambiabilidad: Si bien no es estrictamente parte de SUTVA, ha sido reconocido como uno de los elementos más importantes y estudiados en epidemiología.12 Fuera de ensayos aleatorios, la intercambiabilidad será condicional a los confusores medidos o teorizados. Implica preguntarse: ¿Son comparables los grupos, dada la confusión?, ¿Los tratados y no tratados habrían experimentado en promedio el mismo riesgo de muerte de haber recibido el mismo nivel de tratamiento?, ¿Sigue siendo informativo para el resultado observado saber si son tratados o no, condicional a los niveles de confusión?, ¿Los resultados son independientes de la asignación a tratamiento dado que se tienen las covariables?, ¿Hay independencia entre la asignación al tratamiento y los resultados potenciales?14. Sesgos como el error de medición, variables no-observadas también podrían afectar el cumplimiento de estos supuestos. Sin la intercambiabilidad, el efecto puede ser interpretado como un efecto local, por lo que es un requisito de la validez interna.12
Ajuste por confusión
Incluso si se cumplen los supuestos antes mencionados, es necesario seguir algunas estrategias analíticas para identificar un efecto causal propiamente tal de una exposición/tratamiento (\(A\)) en un resultado \(Y\). En las investigaciones suelen existir vías no-casuales o ‘caminos de puerta trasera’ abiertos que puedan indicar una asociación espuria18. Con todo, la manera que los manejamos los confusores (\(L_{it}\)) puede introducir sesgo.3
Hay aproximaciones condicionales y aquellas basadas en al estandarización.
- Estratificación: Son condicionales al estrato y al modelo. Lo anterior puede hacerse centrando en los niveles de las covariables, por ejemplo en una regresión; o restringiendo el análisis a un nivel de una covariable, lo que se conoce como estratificación propiamente tal.3 Por años, la epidemiología ha estado dominada por la aproximación condicional: controlando por estratificación o regresión.19,20 No obstante, la no-colapsabilidad de efectos no lineales como los odds ratio o hazard ratios implican que efectos condicionales pueden diferir de los marginales. Se fijan los valores de las variables para ajustar, y por tanto, distintas variables a ajustar utilizadas arrojarán otra magnitud o dirección del efecto21.
- Estandarización: Como se mencionó anteriormente, apunta a calcular el resultado contrafactual que habría ocurrido en expuestos de no haberlo estado (asumiendo el cumplimiento de los supuestos nombrados arriba). Originalmente utilizada para ajustar mortalidad (ej., por tramos de edad y sexo), su utilización más reciente apunta a utilizar a la población como objetivo. El contrafactual nunca lo observamos directamente, lo estimamos mediante la pseudo-poblacion en que la exposición es independiente de confusores, permitiendo la estimación de parámetros en toda la población de haber recibido el tratamiento y que no dependen de los valores de las otras covariables.3,19,22
Supuestos en estudios observacionales longitudinales
Cuando se quiere analizar efectos longitudinales, no basta con el cumplimiento de los supuestos planteados arriba. Muchas exposiciones en epidemiología y salud pública son tiempo-variables (valores que toman más de un valor en el tiempo).10 Al trabajar con datos de este tipo, se suman otras fuentes de confusión y sesgo. En tal escenario, la estrategia analítica será de mucha importancia para reducir sesgos.23
Por ejemplo, se generan sesgos de selección en el seguimiento, se introducen confusores tiempo-dependientes, algunos producto del tratamiento anterior (treatment confounder feedback) o ajustes que van en un camino desde la exposición al resultado (adjustment for an intermediate or a descendant of an intermediate). Estos elementos derivan en sobreajuste u estratificación por colliders (al ajustar por un camino naturalmente bloqueado al ser efecto común de una o más variables o descendiente de este).23,24 De tal manera, supuestos como la intercambiabilidad, positividad y consistencia deben ser ampliados a una lógica secuencial. Para ejemplificar, basados en un ejercicio de Hernán25:
Show code
translogo25 <-magick::image_read('./_figs/intercambiabilidad_ej.png')%>%
magick::image_transparent('#59897f') %>%
magick::image_transparent('#68b9ab') %>%
magick::image_transparent('#50736c') %>%
magick::image_transparent('#61a494') %>%
magick::image_transparent('#3f635c') %>%
magick::image_transparent('#609c94') %>%
magick::image_transparent('#60bbad') %>%
magick::image_transparent('#35514e') %>%
magick::image_transparent('#6ebdaf') %>%
magick::image_transparent('#6e9b91') %>%
magick::image_transparent('#659e91') %>%
magick::image_transparent('#73ab9e') %>%
magick::image_transparent('#6bb5a8') %>%
magick::image_transparent('#6aa89b')
translogo25
Figure 1: Diagrama estudio hipotético (Fuente: Hernán, Curso EDX Harvard)
Un set mínimo de confusores a ajustar por \(A_0\) y \(L_1\). Por lo cual, el efecto en tiempo 2 puede interpretarse como \(Y^{a0,a1} \perp \!\!\! \perp A_1=a_1|L_1, A_0=a_0\) : Si hay intercambiabilidad en \(A_1\), condicional a \(L_1\); y considerando que también hay intercambiabilidad en el tiempo 0 (\(Y^{a_0,a_1} \perp \!\!\! \perp A_0\)), concluimos que hay intercambiabilidad secuencial. De igual manera, la consistencia debiese quedar representada como: \[ E(Y|A_0=a_0,A_1=a_1)=E(Y^{a0,a1}|A_0=a_0,A_1=a_1)\], por lo que la respuesta observada según el valor de exposición debiese ser igual al resultado potencial que sería observado bajo el mismo valor de exposición. Por último, la positividad también debiese ser secuencial:16 \[0 < P(A_1 = 1|L_1 = l_1, A_0 = a_0) < 1\,\,\,\,\,\,\&\,\,\,\,\,\,0 < P(A_0 = 1) < 1.\]
“In longitudinal studies with time-dependent confounding, identifying the structure allows us to detect situations in which stratification-based methods would adjust for confounding at the expense of introducing selection bias” (p. 622)23
Marginal structural models
El nombre “marginal structural” apunta a estimar los efectos causales (estructurales) y marginales porque estiman el efecto en toda la población, no sólo condicional a los tratados (average treatment effects o ATE por sus siglas en inglés). Todos ellos pueden estimar el efecto de exposición variable en el tiempo en presencia de confusores tiempo-dependientes afectados por la exposición previa.10,26–31
Se distinguen 3 modelos principales:
- IPTW: El más simple e intuitivo es la ponderación por la probabilidad inversa de asignación a tratamiento (inverse probability treatment weighting o IPTW, por sus siglas en inglés). Es un modelo en dos etapas, en que primero estimamos la probabilidad de recibir tratamiento en base a un vector de confusores \(L\). Luego en una segunda etapa se estima el modelo de regresión como tal, pero ya no existiendo la asociación entre confusor y tratamiento, ya que la pseudopoblación ponderada se encuentra balanceada en estos términos. En otras palabras, la pseudo-población es independiente de la asignación a tratamiento. Sus ventajas radican en lo intuitiva de su utilización, ya que permitiría ver los pesos asignados directamente, siendo además fácil de testear, permite manejar censura longitudinal y puede usarse directamente para MSMs más complejos. Lo problemático de este modelo es que algunas veces los tratados y no tratados tienen características muy distintas en los niveles de L (falta de superposición o lack of common support), lo que puede generar problemas para asignar pesos (ej., asignaría pesos muy extremos e irreales en la práctica). También se ha criticado que es ineficiente en términos computacionales y puede ser inestable además en muestras finitas, además de producir intervalos de confianza más amplios. Para estimar los efectos marginales en muestras correlacionadas, se le utiliza en conjunto con las llamadas ecuaciones de estimación generalizadas (GEE). La distribución marginal de la variable respuesta es modelada como función de las covariables, asumiendo un modelo lineal generalizado o de cuasi-verosimilitud, y la asociación es incorporada para obtener estimaciones más confiables de los parámetros marginales de regresión.10,12
A continuación se presentan los métodos G. Se les denomina métodos G porque pueden ser usados para cualquier tratamiento de manera generalizada, ya sea fijado en el tiempo o variable. Generalizan la primera ecuación (la probabilidad condicional de asignación a tratamiento/exposición) para permitir ajustes por confusores tiempo-dependientes o intermediarios entre exposición temprana y resultado final. Ambos trabajan con un modelo de asignación a tratamiento (pi-model) y un modelo para la variable de resultado (Q-model).
- Fórmula G paramétrica:
\[ \sum_{\bar{l}} E[Y | \bar{A}= \bar{a}, \bar{L} = \bar{l}] ∗ \prod^K_{k= 0}f(l_k|\bar{a}_{k− 1},\bar{l}_{k− 1})\]32
A diferencia de IPTW usa la probabilidad de tratamiento A dadas las covariables L, esta técnica analítica usa la probabilidad de la covariable \(L\) y la probabilidad condicional del resultado \(Y\) dado \(A\) y \(L\). Dicho esto, requiere la esperanza del resultado \(Y\) condicional al tratamiento y confusión en el Tiempo 1 y 2 (exposición e historias de confusión como predictores), más \(Q_1\) equivalente a la esperanza condicional de \(Q_2\) dado \(L1\) (las historias de exposición y confusión como predictor). \[Q_2 = E(Y|L_1,A_1=1, L_2, A_2=1)\] \[Q_1 = E(Q_2|L_1, A_1=1)\]
La densidad conjunta de los datos observados se obtiene derivando por medio de derivaciones Monte Carlo sobre los confusores de \(Q_1\), y luego sustituyéndolos en \(Q_2\) en función del regimen de exposición de interés. Si el modelo de exposición es dinámico se necesita incorporar un modelo de exposición (\(Q_3\)). Los intervalos de confianza se estiman mediante un remuestreo bootstrap paramétrico. El efecto causal marginal se formula de la siguiente manera:
\[\Psi = E(Y^{a_0,a_1}-Y^{0,0})\]
Una de sus desventajas es su sensibilidad a muestras excesivamente grandes en que tiende a rechazar \(H_0\), lo que se ha denominado G-null paradox y razón por la que se recomienda utilizarse sólo en intervenciones en que se sospecha \(H_1\).
- Estimación G:
El primero es un modelo causal que incluye la variable causal de interés y enlaza el resultado contrafactual bajo no-exposición en el seguimiento (ej., el resultado que hubiese ocurrido si no se hubiese observado exposición). La segunda es una regresión logística para predecir la exposición en cada visita basado en la exposición previa y las historias de covariables y el resultado contracfactual. Luego la estimación G busca iterativamente los valores de la variable causal en el primer modelo causal, el que cuando es usado en el segundo modelo, hace la exposición independiente del resultado contrafactual estimado dado los tratamientos previos y la historia de confusores. La estimación G logra ajustar por confusores tiempo-dependientes que son afectados por la exposición previa, ya que examina por separado la asociación entre la exposición y el resultado contrafactual en cada visita y ajustando sólo por los valores de confusores tiempo-dependientes en tiempos anteriores. Los intervalos de confianza pueden construirse en un procedimiento basado en tests o remuestreo Bootstrap no-paramétrico.
Estandariza el riesgo de la distribución conjunta (joint distribution) de los confusores en todo el seguimiento densidad conjunta de los datos observados para generar resultados potenciales bajo distintos escenarios de exposición. Tiene la ventaja de que no asume una distribución paramétrica y por tanto es insensible a una eventual mala especificación.
La gran ventaja de esta herramienta es que puede usarse incluso cuando no hay “common support” (hay un nivel de confusión en que todos los participantes están expuestos o no expuestos y por tanto el supuesto de positividad peligra o near-positivity violations). No obstante y como ya se mencionó, es muy sensible a la mala especificación, lo que llevará a estimaciones sesgadas mientras más etapas (\(k\)) se tenga. Por otra parte, es muy intensiva computacionalmente y requiere mayor programación y datos para resultados estables. Los intervalos de confianza deben obtenerse mediante remuestreos de tipo bootstrap.
Objetivos
Este proyecto apunta a servir como un material introductorio sobre modelos estructurales marginales para el control de sesgos en estudios observacionales con factores de riesgo y/o exposición tiempo-dependientes.33 Para ello se muestra su estimación paso a paso y mediante un ejercicio reproducible. De igual manera, se utiliza una base de datos simulada para ilustrar el proceso.
Ejercicio de simulación
Se generó una base de datos de 1.000 observaciones con las siguientes características en base al ejemplo de Bounthavong,34 con la finalidad de estimar el efecot de tratamiento residencial en el recuento de meses libres de readmisión:
t1= Tratamiento residencial Tiempo 0 ( \(A_0\) ) ∼ Bernoulli(p=.23)
t2= Tratamiento residencial Tiempo 1 ( \(A_1\) ) ∼ Bernoulli(p**=.22)
t3= Tratamiento residencial Tiempo 2 ( \(A_2\) ) ∼ Bernoulli(p**=.15)
c1= Policonsumo ∼ Bernoulli(p**=.32)
c2 = Edad de inicio consumo de sustancias [
edad_ini] ∼ N(3,0.5)dt1 = Meses sin readmisión Primer tto. ∼ Poisson(λ**=4)
dt2 = Meses sin readmisión Segundo tto. ∼ Poisson(λ**=3)
dt3 = Meses sin readmisión Tercer tto. ∼ Poisson(λ**=2)
v1 = Confusor tiempo-dependiente (t1) ∼ \(5+ 0.4\beta_{Tto\,residencial\,basal} +0.78\beta_{Meses\,sin\,readmision\,1}+ N(0, \sqrt{0.99})\)
v2 = Confusor tiempo-dependiente (t2) ∼ \(5+ 0.4\beta_{Tto\,residencial\,(2^{do})} +0.78\beta_{Meses\,sin\,readmision\,1}+ N(0, \sqrt{0.55})\)
v3 = Confusor tiempo-dependiente (t3) ∼ \(5+ 0.4\beta_{Tto\,residencial\,(3^{er})} +0.78\beta_{Meses\,sin\,readmision\,1}+ N(0, \sqrt{0.33})\)
Lo anterior, aplicando un modelo de ponderación inversa de la
probabilidad de asignación de tratamiento, y luego un modelo marginal
estructural en base a un modelo – GEE (generalized estimating ecuations)
con errores estándar robustos poisson, con estructura de correlación
autorregresiva de orden 1 (de tal manera que \(E[Tratamiento=1, tiempo=1 | Tratamiento=1,
tiempo=(tiempo-1)]\)) para datos anidados por paciente (que
permita reflejar la dependencia en los resultados del tratamiento
anterior). Se comparan los resultados de un modelo que no consideró la
ponderación, versus uno que sí consideró la ponderación.
Los datos poseen la siguiente estructura:
Show code
#Semilla para replicar los resultados
set.seed(2125)
#definir confusor c1 (policonsumo al ingreso==1)
policonsumo <-rbinom(1000,1,0.32)
#definir confusor c2 (edad de inicio de consumo)
edad_ini <-round(exp(rnorm(1000, 3, 0.5)))
#definir tratamiento residencial al tiempo 0
t1 <-rbinom(1000,1,0.23)
#definir tratamiento residencial al tiempo 1
t2 <-rbinom(1000,1,0.22)
#definir tratamiento residencial al tiempo 2
t3 <-rbinom(1000,1,0.15)
#meses libre de readmisión T1
dt1 <-rpois(1000,4)
#meses libre de readmisión T2
dt2 <-rpois(1000,3)
#meses libre de readmisión T3
dt3 <-rpois(1000,1)
#definir el confusor tiempo-dependiente v1 como una función de t1 y dt1
v1 <- (0.4*t1 + 0.78*dt1 + rnorm(1000, 0, sqrt(0.99))) + 5
#definir el confusor tiempo-dependiente v2 como una función de t2 y dt2
v2 <- (0.4*t2 + 0.78*dt2 + rnorm(1e3, 0, sqrt(0.55))) + 5
#definir el confusor tiempo-dependiente v3 como una función de t3 y dt3
v3 <- (0.4*t3 + 0.78*dt3 + rnorm(1e3, 0, sqrt(0.33))) + 5
datos <- data.frame(policonsumo, edad_ini, v1, v2, v3, t1, t2, t3, dt1, dt2, dt3) %>%
# Generamos una columna id para cada sujeto
dplyr::mutate(id = row_number()) %>%
# Convertimos la base de datos desde formato ancho a largo de todas las variables tiempo-dependientes. Aclaramos que la estructura de los nombres lleva inscrito primero el valor de interés, seguido por el tiempo/ola de aplicación. Y obtenemos el tiempo tomando el valor dígito de cada columna
tidyr::pivot_longer(cols=c("v1", "v2", "v3", "t1", "t2", "t3", "dt1", "dt2", "dt3"), names_to = c('.value',"time"), names_pattern = "(.*)(\\d+)$") %>%
dplyr::rename("c1"="policonsumo","c2"="edad_ini") %>%
# Ordenamos la base de datos por id y tiempo/ola
dplyr::arrange(id, time) %>%
# reordenamos las columnas
dplyr::select(id, time, c1, c2, t, v, dt) %>%
# formateamos la base de datos como data.frame
as.data.frame() %>%
dplyr::mutate(time=as.numeric(time))
#Primeros 4 filas
head(datos,4) %>%
knitr::kable("html", caption="Muestra de la base de datos (Primeras 4 filas)") %>%
kableExtra::kable_classic(bootstrap_options = c("striped", "hover"),font_size = 10)
| id | time | c1 | c2 | t | v | dt |
|---|---|---|---|---|---|---|
| 1 | 1 | 0 | 30 | 0 | 10.633698 | 7 |
| 1 | 2 | 0 | 30 | 0 | 6.549493 | 3 |
| 1 | 3 | 0 | 30 | 0 | 7.595415 | 3 |
| 2 | 1 | 1 | 44 | 0 | 9.276188 | 2 |
Se genera un diagrama causal para reflejar la estructura de confusión y sesgo presente en los datos en caso de estimar el efecto de \(X\) en \(Y\).
Show code
datos_mod<-
#transformamos la base de datos a formato ancho, en que cada fila es un sujeto
pivot_wider(dplyr::mutate(datos,time=time-1), id_cols=c("id","c1","c2"),names_from="time", values_from=c("t","v","dt"))%>%
dplyr::mutate(across(paste0("v_",0:2),~as.integer(.x*100))) %>%
dplyr::mutate(c2=as.integer(c2)) %>%
dplyr::rename("o_0"="c1","p_0"="c2") %>%
dplyr::mutate(o_1=o_0,o_2=o_1,p_1=p_0,p_2=p_1) %>%
dplyr::relocate(o_0, .before=o_1) %>%
dplyr::relocate(p_0, .before=p_1) %>%
data.frame()
#Agregamos la historia de la exposición
datos_mod2<-makehistory.one(datos_mod,id="id",exposure="t",times=0:2)
d1_noc<-
diagnose(input=datos_mod2,
diagnostic=1,
censoring="no",
approach="none",
scope="all",
id="id",
times.exposure=c(0,1,2),
times.covariate=c(0,1,2),
exposure="t",
temporal.covariate=c("v"),
static.covariate=c("o","p"),
history="h",
ignore.missing.metric="no",
loop="yes",
sd.ref="no")
# d1_noc %>% ggplot(aes(x = name.cov, y = abs(SMD))) +
# geom_hline(yintercept = 0.2, color = "black", size = 0.2, linetype="dashed") +
# coord_flip() +
# theme_classic() +
# theme(legend.key = element_blank()) +
# labs(y = "Diferencias promedio estandarizadas en medias",x = "",color="",shape="") +
# geom_point()+
# facet_wrap(time.covariate+time.exposure~H, ncol=6)+
# ggtitle("Tiempo de exposición (vertical), historia de exposición (horizontal)")+
# theme(panel.spacing.x = unit(0,"line"))
d1_noc %>% knitr::kable("html", caption="Balance covariables confusoras", col.names = c("E","Historia de exposición","Covariable","Tiempo exposición","Tiempo covariable","D", "SMD", "N", "N° Expuestos")) %>%
kableExtra::kable_classic() %>%
kableExtra::scroll_box(width = "100%", height = "250px")
| E | Historia de exposición | Covariable | Tiempo exposición | Tiempo covariable | D | SMD | N | N° Expuestos |
|---|---|---|---|---|---|---|---|---|
| 1 | H | o | 0 | 0 | -0.0004677 | -0.0010017 | 1000 | 219 |
| 1 | H | p | 0 | 0 | -0.2377060 | -0.0205881 | 1000 | 219 |
| 1 | H | v | 0 | 0 | 38.1333088 | 0.2079093 | 1000 | 219 |
| 1 | H0 | o | 1 | 0 | -0.0376727 | -0.0806969 | 781 | 182 |
| 1 | H0 | p | 1 | 0 | 1.2033976 | 0.1038620 | 781 | 182 |
| 1 | H0 | v | 1 | 0 | -0.4607037 | -0.0025040 | 781 | 182 |
| 1 | H0 | v | 1 | 1 | 38.8674439 | 0.2543879 | 781 | 182 |
| 1 | H00 | o | 2 | 0 | 0.0640780 | 0.1363079 | 599 | 78 |
| 1 | H00 | p | 2 | 0 | -1.3770117 | -0.1218171 | 599 | 78 |
| 1 | H00 | v | 2 | 0 | 28.3163049 | 0.1520682 | 599 | 78 |
| 1 | H00 | v | 2 | 1 | -29.7416950 | -0.1921208 | 599 | 78 |
| 1 | H00 | v | 2 | 2 | 34.2484620 | 0.3483719 | 599 | 78 |
| 1 | H01 | o | 2 | 0 | -0.0705128 | -0.1545784 | 182 | 26 |
| 1 | H01 | p | 2 | 0 | -0.3333333 | -0.0266485 | 182 | 26 |
| 1 | H01 | v | 2 | 0 | 28.9423077 | 0.1641306 | 182 | 26 |
| 1 | H01 | v | 2 | 1 | -28.2756410 | -0.1946057 | 182 | 26 |
| 1 | H01 | v | 2 | 2 | 12.1602564 | 0.1200610 | 182 | 26 |
| 1 | H1 | o | 1 | 0 | 0.0081679 | 0.0174353 | 219 | 46 |
| 1 | H1 | p | 1 | 0 | -2.5207339 | -0.2218004 | 219 | 46 |
| 1 | H1 | v | 1 | 0 | -59.9141744 | -0.3318739 | 219 | 46 |
| 1 | H1 | v | 1 | 1 | 15.2354863 | 0.0990654 | 219 | 46 |
| 1 | H10 | o | 2 | 0 | 0.0928105 | 0.1985543 | 173 | 20 |
| 1 | H10 | p | 2 | 0 | 8.7666667 | 0.7602868 | 173 | 20 |
| 1 | H10 | v | 2 | 0 | 83.3794118 | 0.4653904 | 173 | 20 |
| 1 | H10 | v | 2 | 1 | -18.5921569 | -0.1234705 | 173 | 20 |
| 1 | H10 | v | 2 | 2 | 21.0398693 | 0.2208975 | 173 | 20 |
| 1 | H11 | o | 2 | 0 | -0.1291291 | -0.2710264 | 46 | 9 |
| 1 | H11 | p | 2 | 0 | -3.0990991 | -0.3307462 | 46 | 9 |
| 1 | H11 | v | 2 | 0 | 69.5165165 | 0.3861858 | 46 | 9 |
| 1 | H11 | v | 2 | 1 | -77.6756757 | -0.4692929 | 46 | 9 |
| 1 | H11 | v | 2 | 2 | 61.5975976 | 0.6501729 | 46 | 9 |
A partir de la tabla anterior, se muestran algunos desbalances entre las covariables según historia de exposición, tiempo de exposición y covariables. Estos desbalances exigen ajustar por algunas variables.35
Show code
#```{r, echo=F, dev.args = list(bg = 'transparent'), fig.align="center", out.width=550, error=T}
#fig.width = 7, fig.height = 5 muy grande
#out.height=450,
dag34 <- dagitty('dag {
bb="0,0,1,1"
L1_c1 [pos="0.117,0.205"]
L2_c2 [pos="0.222,0.112"]
dt1 [outcome,pos="0.294,0.491"]
dt2 [outcome,pos="0.572,0.491"]
dt3 [outcome,pos="0.838,0.491"]
t1 [exposure,pos="0.092,0.491"]
t2 [exposure,pos="0.380,0.491"]
t3 [exposure,pos="0.657,0.491"]
v1 [pos="0.289,0.686"]
v2 [pos="0.576,0.668"]
v3 [pos="0.849,0.657"]
L1_c1 -> dt1
L1_c1 -> dt2
L1_c1 -> dt3
L1_c1 -> t1
L1_c1 -> t2
L1_c1 -> t3
L2_c2 -> dt1
L2_c2 -> dt2
L2_c2 -> dt3
L2_c2 -> t1
L2_c2 -> t2
L2_c2 -> t3
t1 -> dt1
t1 -> t2 [pos="0.210,0.652"]
t2 -> dt2
t2 -> t3 [pos="0.514,0.677"]
t3 -> dt3
v1 -> dt1
v1 -> t1
v1 -> t2
v2 -> dt2
v2 -> t2
v2 -> t3
v3 -> dt3
v3 -> t3
dt1 -> dt2
dt2 -> dt3
}')
tidy_dag34 <- tidy_dagitty(dag34) %>%
dplyr::mutate(label=dplyr::case_when(grepl("A0",as.character(name))~"A0",
T~as.character(name))) %>%
dplyr::mutate(label2=dplyr::case_when(grepl("L|v",name)~"adj",grepl("LM",name)~"white",
T~"black")) %>%
dplyr::mutate(adjusted=factor(dplyr::case_when(grepl("L|v",name)~"adjusted",T~"unadjusted")))
edge_function <- ggdag:::edge_type_switch("link_arc")
dag34_plot<-
ggdag:::if_not_tidy_daggity(tidy_dag34) %>% ggdag:::node_status() %>%
ggplot2::ggplot(ggplot2::aes(x = x, y = y, xend = xend,
yend = yend, color = status, shape=factor(adjusted)))+
#edge_function()+
scale_adjusted()+ ggdag:::breaks(c("exposure", "outcome","latent"))+
theme_dag()+
geom_dag_edges_arc(curvature = c(rep(0,12),5,5,0,5,0,5,rep(0,9)))+ #14 y 16 de 24
ggdag:::geom_dag_point(size = 16)+
ggdag:::geom_dag_label_repel(ggplot2::aes_string(label = "label",
fill = "status"), size = 4.88, col = "white",
show.legend = FALSE)+
scale_shape_manual(values = c(15, 16), name="Ajustado", labels=c("Sí", "No"))+
scale_fill_manual(values = c("gray80", "gray30","gray30"), name="Estatus",na.value="black", labels=c("Exposición", "Resultado"), limits = c('exposure', 'outcome'))+
scale_color_manual(values = c("gray75", "gray35","gray30"), name="Estatus",na.value="black", labels=c("Exposición", "Resultado"), limits = c('exposure', 'outcome'))+#E6E6E6
guides(linetype="none", edge_alpha="none", shape="none")+
guides(color=guide_legend(override.aes = list(arrow = NULL)))+#,guide_colourbar(order = 1)
theme(plot.caption = element_text(hjust = 0))+
theme(legend.position = "bottom", aspect.ratio=6/10)+
labs(caption="Nota. Ak= Modalidad (Residencial/Ambulatoria); dtk= Meses libre de readmisión;\nvk= Confusor tiempo-dependiente; c1= Policonsumo; c2= Edad de inicio consumo de sustancias")
dag34_plot
Figure 2: Diagrama de la estructura causal del ejercicio de simulación
Para simplificar el ejemplo, se asume que todos cursaron el tratamiento y no hay censura, que la velocidad es constante y los riesgos son constantes y proporcionales.36 a diferencia de los modelos Cox en que el parámetro lambda (tasa de intensidad) sigue un proceso estocástico.37
Los ponderadores de asignación a tratamiento aquí utilizada tiene en cuenta una estructura de exposición y confusión tiempo-dependientes. Para cada unidad observada, esta función genera una ponderación de probabilidad inversa de tratamiento para cada etapa de seguimiento. Para ello se incluye en el numerador los confusores invariantes y el denominador todos los confusores de interés, incluidos los tiempo-dependientes.38
Show code
#Paquetes necesarios para generar los análisis
library(geepack)
library(survey)
library(ipw)
library(reshape)
library(MuMIn)
#Estimamos el modelo de probabilidad inversa con el numerador (elementos invariantes en el tiempo) y el denominador todos los confusores de interés.
w <- ipwtm(
exposure = t,
family = "binomial",
link = "logit",
# Invariantes en el tiempo
numerator = ~ factor(c1) + c2,
# Todos los confusores
denominator = ~ v + factor(c1) + c2,
id = id,
timevar=time,
type="first",
data = datos)
summary(w$ipw.weights)
#Se incorporan los pesos a la base de datos
datos<-
dplyr::bind_cols(datos, ipw=w$ipw.weights)
###########################################################################################
# Para obtener las salidas de los coeficientes
###########################################################################################
summary.geem <- function(object, ...) {
Coefs <- matrix(NA,nrow=length(object$beta),ncol=5)
Coefs[,1] <- c(object$beta)
naive <- is.character(object$var)
if(!naive && any(diag(object$var) < 0) ){
naive <- TRUE
warning("Some elements of robust variance estimate < 0. Reporting model based SE.")
}
Coefs[,2] <- sqrt(diag(object$naiv.var))
if(naive){Coefs[,3] <- rep(NA, length(object$beta))}else{Coefs[,3] <- sqrt(diag(object$var))}
if(naive){Coefs[,4] <- Coefs[,1]/Coefs[,2]}else{Coefs[,4] <- Coefs[,1]/Coefs[,3]}
Coefs[,5] <- round(2*pnorm(abs(Coefs[,4]), lower.tail=F), digits=8)
colnames(Coefs) <- c("Estimates","Model SE","Robust SE", "wald", "p")
summ <<- list(beta = Coefs[,1], se.model = Coefs[,2], se.robust = Coefs[,3], wald.test = Coefs[,4], p = Coefs[,5],
alpha = object$alpha, corr = object$corr, phi = object$phi, niter = object$niter, clusz = object$clusz,
coefnames = object$coefnames, weights=object$weights, biggest.R.alpha = object$biggest.R.alpha)
class(summ) <- 'summary.geem'
return(summ)
}
###########################################################################################
# Primer modelo GEE - GEE con errores estándar robustos (no IPTW)
###########################################################################################
mod1<-geem(dt~ t+ time+ factor(c1)+ c2, id=id ,data = data.frame(datos), family=poisson,
corstr="ar1")
summary(mod1)
summ_mod1<-
cbind.data.frame(cc=summ$coefnames,
Estimate=summ$beta,Std.err=summ$se.robust) %>%
dplyr::mutate(lwr=Estimate-qnorm((1+0.95)/2)*Std.err,upr=Estimate+qnorm((1+0.95)/2)*Std.err) %>%
dplyr::mutate(across(where(is.numeric),~round(exp(.),2)))
invisible("Quasipoisson")
Si bien se probaron distribuciones binomiales negativas para ver el grado de sobredispersión de los datos y se comparan mediante criterios de parsimonia bajo cuasi-verosimilitud (QIC), de antemano sabemos que no habrá ni sobre ni subdispersión dado que son datos simulados.
Show code
#, results='hide'
#https://gist.github.com/avallecam/56af06f46e5544c3af0f46344df20989
#https://stackoverflow.com/questions/49089476/any-updates-to-model-negative-binomial-distribution-data-with-gee-in-r
#https://www.rdocumentation.org/packages/geepack/versions/1.3.4/topics/geese
#https://stackoverflow.com/questions/13946540/negative-binomial-in-gee?rq=1
#_#_#_#_ POIS
mod1b<- geem(dt~ t+ time+ factor(c1)+ c2, id=id ,weights=ipw, data = data.frame(datos), family=poisson,
corstr="ar1")
invisible(summary.geem(mod1b)) #para ver los resultados
summ_mod1b<-
cbind.data.frame(cc=summ$coefnames,
Estimate=summ$beta,Std.err=summ$se.robust) %>%
dplyr::mutate(lwr=Estimate-qnorm((1+0.95)/2)*Std.err,upr=Estimate+qnorm((1+0.95)/2)*Std.err) %>%
dplyr::mutate(across(where(is.numeric),~round(exp(.),2)))
#_#_#_#_ NB2
mod2b<- geem(dt~ t+ time+ factor(c1)+ c2, id=id, weights=ipw, data = datos,
family = MASS::negative.binomial(4^2), corstr = "ar1") #El valor conocido del parámetro theta
invisible(summary.geem(mod2b)) #para ver los resultados
summ_mod2b<-
cbind.data.frame(cc=summ$coefnames,
Estimate=summ$beta,Std.err=summ$se.robust) %>%
dplyr::mutate(lwr=Estimate-qnorm((1+0.95)/2)*Std.err,upr=Estimate+qnorm((1+0.95)/2)*Std.err) %>%
dplyr::mutate(across(where(is.numeric),~round(exp(.),2)))
#_#_#_#_ NB1
mod3b<- geem(dt~ t+ time+ factor(c1)+ c2, id=id, weights=ipw, data = datos,
family = MASS::negative.binomial(3^2), corstr = "ar1")
invisible(summary.geem(mod3b)) #para ver los resultados
summ_mod3b<-
cbind.data.frame(cc=summ$coefnames,
Estimate=summ$beta,Std.err=summ$se.robust) %>%
dplyr::mutate(lwr=Estimate-qnorm((1+0.95)/2)*Std.err,upr=Estimate+qnorm((1+0.95)/2)*Std.err) %>%
dplyr::mutate(across(where(is.numeric),~round(exp(.),2)))
#_#_#_#_ NB 0.5
mod4b<- geem(dt~ t+ time+ factor(c1)+ c2, id=id, weights=ipw, data = datos,
family = MASS::negative.binomial(2^2), corstr = "ar1")
invisible(summary.geem(mod4b)) #para ver los resultados
summ_mod4b<-
cbind.data.frame(cc=summ$coefnames,
Estimate=summ$beta,Std.err=summ$se.robust) %>%
dplyr::mutate(lwr=Estimate-qnorm((1+0.95)/2)*Std.err,upr=Estimate+qnorm((1+0.95)/2)*Std.err) %>%
dplyr::mutate(across(where(is.numeric),~round(exp(.),2)))
#Poisson tiene una varianza igual al promedio, por lo que el valor theta es el denominador de la media al cuadrado y da una impresión del grado de sobredispersión que hay en los datos
#Se muestra el criterio de cuasi-verosimilitud bajo el modelo de independencia (QIC; Pan, 2001) para ver el grado de ajuste del modelo.
model.sel(mod1b, mod2b, mod3b,mod4b, rank = QIC)
Model selection table
(Int) c2 fct(c1) t tim family qLik
mod1b 2.139 -0.0005898 + -0.1268 -0.6095 p(l) 771.262
mod2b 2.157 -0.0005834 + -0.1263 -0.6200 NB(16,l) 771.177
mod3b 2.169 -0.0005770 + -0.1260 -0.6266 NB(9,l) 752.003
mod4b 2.195 -0.0005573 + -0.1253 -0.6409 NB(4,l) 659.108
QIC delta weight
mod1b -1536.6 0.00 0.561
mod2b -1536.1 0.49 0.439
mod3b -1497.5 39.06 0.000
mod4b -1311.2 225.37 0.000
Abbreviations:
family: NB(16,l) = 'Negative Binomial(16,log)',
NB(4,l) = 'Negative Binomial(4,log)',
NB(9,l) = 'Negative Binomial(9,log)', p(l) = 'poisson(log)'
Models ranked by QIC(x) Como era de esperarse, el modelo mod1b asumiendo una
distribución Poisson tiene un mejor ajuste que el resto (asumiendo
distribución binomial negativa con valor \(\theta\) \(4^2\), \(3^2\) y \(2^2\), respectivamente). Una vez que se
genera el modelo GEE sin ajustar por la probabilidad inversa de
asignación a tratamiento, se contrasta con el modelo en el que se
pondera por dicha probabilidad.
Show code
mt<-
rbind.data.frame(
summ_mod1[which(summ_mod1$cc=="t"),],
summ_mod1b[which(summ_mod1b$cc=="t"),]
)
rownames(mt) <- c("Modelo sin IPTW", "Modelo con IPTW")
mt %>%
dplyr::select(-cc,-`Std.err`) %>%
dplyr::mutate(paste0(lwr,", ",upr)) %>%
dplyr::select(-lwr,-upr) %>%
knitr::kable("html", caption="Comparación resultados", col.names=c("RR", "95% IC")) %>%
kableExtra::kable_classic(bootstrap_options = c("striped", "hover"),font_size = 11) %>%
kableExtra::add_footnote("Nota. Intervalo de confianza al 95%")
| RR | 95% IC | |
|---|---|---|
| Modelo sin IPTW | 0.99 | 0.94, 1.05 |
| Modelo con IPTW | 0.88 | 0.83, 0.94 |
| a Nota. Intervalo de confianza al 95% |
El modelo generalizado sin ajustar señala que no hay diferencias en el número de meses libre de reingresar a tratamiento entre los tratamientos residenciales y ambulatorios, ya que los intervalos de confianza se superponen al valor nulo (RR= 0.99 95% IC: 0.94, 1.05). No obstante, el modelo que presenta ponderación por el tratamiento indica que los tratamientos residenciales presentan 1/0.88 o un riesgo 14% menor de estar más meses libre de readmisión, por lo que serían readmitidos 14% más rápido (RR= 0.88 95% IC: 0.83, 0.94). Dicho lo anterior, estar en tratamiento residencial tendría un efecto perjudicial una vez ponderando por la probabilidad inversa de asignación a tratamiento en el modelo estructural marginal.
Discusión
Este estudio fue una aproximación muy preliminar al uso de modelos marginales estructurales. En el ejemplo presentado con datos simulados, se pudo apreciar cómo el uso de modelos marginales reducen la confusión generada por variables tiempo-dependientes y reproduce de mejor manera el proceso de generación de datos especificado en la simulación.
Es necesario tener en consideración que hay una tensión entre definir tan fielmente el tratamiento para cumplir SUTVA y tener la capacidad de generalizar un constructo más allá de la operacionalización de una exposición, por lo que deben tenerse en cuenta tales consideraciones al momento de diseñar estudios observacionales con objetivos de inferencia causal.12
Por otra parte, es aconsejable que al momento de elegir la herramienta analítica se tenga en cuenta sus ventajas y desventajas. En la práctica, los modelos g-paramétricos y la estimación g pueden ser difíciles de implementar, particularmente cuando hay muchos tiempos de observación o cuando el tratamiento y los confusores son tiempo-dependientes en varias etapas.
Con todo, al día de hoy se están produciendo grandes esfuerzos para superar la barrera computacional que dificulta la aplicación de estos modelos.39 Y si por el momento no hay versiones del todo consolidadas de métodos-g para observar tratamientos con tiempo continuo y resultados de supervivencia censurados (más allá de IPTW más la probabilidad inversa de censura), de a poco se han ido publicando algunos estudios prometedores en esta materia.40
Limitaciones
En primer lugar, existen los modelos estructurales anidados (Structural Nested Models), que apuntan a identificar el efecto condicional a un regímenes de tratamiento específicos.28 Estos modelos no fueron abordados en este estudio.
Por otra parte, particularmente en datos en que la variable de interés mide una duración o una supervivencia hasta que se produce un evento, pudiesen existir confusores tiempo-modificados en que no sólo adoptan distintos valores en el tiempo, si no que su efecto va cambiando conforme pasa el tiempo (time-modified confounding).29 Dada su complejidad, dicha aproximación va más allá del alcance de este estudio.
De igual manera, existen modelos doble o triplemente robustos que permiten combinar la probabilidad inversa de asignación a tratamiento con métodos g.19
Por último, se han publicado varias herramientas de aprendizaje automático o Machine Learning (bagging, boosting, random forests, neural networks, SuperLearner), que son más capaces de generar un espacio de comparación razonable entre tratados/expuestos y no tratados/no expuestos (no asigna ponderadores muy extremos, penalizan por escasez de superposición en los niveles de una covariable, entre otras ventajas). Con todo y dada la extensión de este trabajo, no han podido ser abordados en el presente estudio.
Información de la sesión
Show code
Sys.getenv("R_LIBS_USER")
[1] "H:/Mi unidad/PERSONAL ANDRES/UCH_salud_publica/asignaturas/0_TESIS/renv/library/R-4.1/x86_64-w64-mingw32;C:/Program Files/R/R-4.1.2/library"Show code
R version 4.1.2 (2021-11-01)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)
Matrix products: default
locale:
[1] LC_COLLATE=Spanish_Chile.1252 LC_CTYPE=Spanish_Chile.1252
[3] LC_MONETARY=Spanish_Chile.1252 LC_NUMERIC=C
[5] LC_TIME=Spanish_Chile.1252
attached base packages:
[1] grid stats graphics grDevices datasets utils
[7] methods base
other attached packages:
[1] confoundr_1.2 MuMIn_1.46.0
[3] geeM_0.10.1 reshape_0.8.8
[5] ipw_1.0-11 survey_4.1-1
[7] Matrix_1.4-1 geepack_1.3.4
[9] ggdag_0.2.4 dagitty_0.3-1
[11] qrcode_0.1.4 distill_1.4
[13] ggiraph_0.7.10 choroplethrMaps_1.0.1
[15] choroplethr_3.7.0 acs_2.1.4
[17] XML_3.99-0.9 choroplethrAdmin1_1.1.1
[19] chilemapas_0.2 sf_1.0-7
[21] compareGroups_4.5.1 mstate_0.3.2
[23] pdftools_3.1.1 flexsurv_2.2
[25] timereg_2.0.1 future_1.23.0
[27] coxphw_4.0.2 survival_3.2-11
[29] cowplot_1.1.1 psych_2.1.9
[31] forcats_0.5.1 stringr_1.4.0
[33] dplyr_1.0.7 purrr_0.3.4
[35] readr_2.1.0 tidyr_1.1.4
[37] tibble_3.1.6 ggplot2_3.3.5
[39] tidyverse_1.3.1 datapasta_3.1.0
[41] widgetframe_0.3.1 htmlwidgets_1.5.4
[43] fontawesome_0.2.2 xaringanthemer_0.4.1
[45] xaringanExtra_0.5.5 xaringan_0.23
[47] DiagrammeR_1.0.9 bibtex_0.4.2.3
[49] magick_2.7.3 pagedown_0.17.1
[51] RefManageR_1.3.0 ggpattern_0.4.2
[53] showtext_0.9-5 showtextdb_3.0
[55] sysfonts_0.8.8 here_1.0.1
[57] devtools_2.4.2 usethis_2.1.3
[59] pacman_0.5.1
loaded via a namespace (and not attached):
[1] rappdirs_0.3.3 downlit_0.4.0 muhaz_1.2.6.4
[4] R.methodsS3_1.8.2 knitr_1.36 R.utils_2.12.0
[7] data.table_1.14.2 rpart_4.1-15 generics_0.1.3
[10] callr_3.7.0 mice_3.14.0 proxy_0.4-26
[13] chron_2.3-56 tzdb_0.2.0 webshot_0.5.2
[16] xml2_1.3.2 lubridate_1.8.0 assertthat_0.2.1
[19] viridis_0.6.2 xfun_0.31 hms_1.1.1
[22] jquerylib_0.1.4 evaluate_0.14 fansi_0.5.0
[25] dbplyr_2.1.1 readxl_1.3.1 igraph_1.2.11
[28] DBI_1.1.1 tmvnsim_1.0-2 Rsolnp_1.16
[31] jsonvalidate_1.3.2 stats4_4.1.2 Rmpfr_0.8-9
[34] ellipsis_0.3.2 backports_1.4.0 V8_4.2.0
[37] vctrs_0.3.8 remotes_2.4.1 cachem_1.0.6
[40] withr_2.4.2 ggforce_0.3.3 tigris_1.6.1
[43] HardyWeinberg_1.7.4 checkmate_2.0.0 rgdal_1.5-32
[46] ggmap_3.0.0 prettyunits_1.1.1 mnormt_2.0.2
[49] svglite_2.0.0 cluster_2.1.2 crayon_1.4.2
[52] tidycensus_1.2.2 crul_1.2.0 labeling_0.4.2
[55] pkgconfig_2.0.3 units_0.8-0 tweenr_1.0.2
[58] nlme_3.1-152 pkgload_1.2.3 nnet_7.3-16
[61] rlang_1.0.4 globals_0.14.0 RJSONIO_1.3-1.6
[64] lifecycle_1.0.1 httpcode_0.3.0 modelr_0.1.8
[67] polyclip_1.10-0 cellranger_1.1.0 rprojroot_2.0.2
[70] flextable_0.6.10 boot_1.3-28 reprex_2.0.1
[73] base64enc_0.1-3 processx_3.5.2 png_0.1-7
[76] viridisLite_0.4.0 rjson_0.2.20 bitops_1.0-7
[79] R.oo_1.25.0 KernSmooth_2.23-20 visNetwork_2.1.0
[82] classInt_0.4-3 maptools_1.1-2 parallelly_1.29.0
[85] qpdf_1.1 jpeg_0.1-9 scales_1.1.1
[88] memoise_2.0.1 magrittr_2.0.1 plyr_1.8.6
[91] compiler_4.1.2 kableExtra_1.3.4 RColorBrewer_1.1-2
[94] cli_3.1.0 listenv_0.8.0 ps_1.6.0
[97] htmlTable_2.3.0 Formula_1.2-4 MASS_7.3-54
[100] tidyselect_1.1.2 stringi_1.7.5 highr_0.9
[103] mitools_2.4 yaml_2.2.1 askpass_1.1
[106] ggrepel_0.9.1 latticeExtra_0.6-29 sass_0.4.0
[109] tools_4.1.2 rmapshaper_0.4.6 future.apply_1.8.1
[112] parallel_4.1.2 RgoogleMaps_1.4.5.3 rstudioapi_0.13
[115] uuid_1.0-3 foreign_0.8-81 gridExtra_2.3
[118] farver_2.1.0 ggraph_2.0.5 digest_0.6.28
[121] lava_1.6.10 geojsonlint_0.4.0 quadprog_1.5-8
[124] Rcpp_1.0.7 broom_1.0.0 writexl_1.4.0
[127] httr_1.4.2 gdtools_0.2.3 WDI_2.7.6
[130] colorspace_2.0-2 rvest_1.0.2 fs_1.5.0
[133] truncnorm_1.0-8 splines_4.1.2 statmod_1.4.36
[136] graphlayouts_0.8.0 sp_1.4-6 renv_0.15.5
[139] sessioninfo_1.2.1 systemfonts_1.0.3 gmp_0.6-5
[142] jsonlite_1.7.2 tidygraph_1.2.1 testthat_3.1.0
[145] R6_2.5.1 Hmisc_4.6-0 pillar_1.7.0
[148] htmltools_0.5.2 glue_1.5.0 fastmap_1.1.0
[151] deSolve_1.30 class_7.3-19 codetools_0.2-18
[154] pkgbuild_1.2.0 mvtnorm_1.1-3 utf8_1.2.2
[157] lattice_0.20-44 bslib_0.3.1 numDeriv_2016.8-1.1
[160] curl_4.3.2 officer_0.4.1 zip_2.2.0
[163] rmarkdown_2.14 desc_1.4.0 munsell_0.5.0
[166] e1071_1.7-9 haven_2.4.3 gtable_0.3.0 Warning: H:/Mi unidad/PERSONAL ANDRES/UCH_salud_publica/asignaturas/
0_TESISShow code
save.image("__analisis.RData")
unlink("*_cache", recursive = T, force = T, expand = TRUE)
sesion_info <- devtools::session_info()
dplyr::select(
tibble::as_tibble(sesion_info$packages),
c(package, loadedversion, source)
) %>%
DT::datatable(filter = 'top', colnames = c('Row number' =1,'Variable' = 2, 'Percentage'= 3),
caption = htmltools::tags$caption(
style = 'caption-side: top; text-align: left;',
'', htmltools::em('Paquetes estadísticos utilizados')),
options=list(
initComplete = htmlwidgets::JS(
"function(settings, json) {",
"$(this.api().tables().body()).css({'font-size': '80%'});",
"}")))