example of application in time series analysis - see vignette(“dlnmTS”)
para la aplicación de modelos lineales y no lineales de retardo distribuido (DLM y DLNM) en el análisis de series temporales. El desarrollo de DLM y DLNM y la implementación del software original para datos de series de tiempo se ilustran enGasparrini et al.[2010] yGasparrini[2011].
Los ejemplos descritos en las siguientes secciones cubren la mayoría de las aplicaciones estándar de la metodología DLNM para datos de series de tiempo y exploran las capacidades del paquete dlnm para especificar, resumir y trazar esta clase de modelos. A pesar de la aplicación específica sobre los efectos en la salud de la contaminación del aire y la temperatura, estos ejemplos se generalizan fácilmente a diferentes temas y forman una base para el análisis de este conjunto de datos u otras fuentes de datos de series temporales. Los resultados incluidos en este documento no pretenden representar hallazgos científicos, sino que se informan con el único propósito de ilustrar las capacidades deldlnmpaquete.
2 Datos
Los ejemplos incluidos en la las asociaciones entre la contaminación del aire y la temperatura con la mortalidad, utilizando un conjunto de datos de series temporales con observaciones diarias para la ciudad de Chicago en el período 1987-2000.
Este conjunto de datos está incluido en el paquete como el marco de datoschicago nmmaps,y se describe en la página de ayuda relacionada (help(chicagoNMMAPS)y la viñetadlnmResumen). Después de cargar el paquete en elRsesión, echemos un vistazo a las tres primeras observaciones:
library(dlnm)## This is dlnm 2.4.7. For details: help(dlnm) and vignette('dlnmOverview').library(DT)
DT::datatable(chicagoNMMAPS)El conjunto de datos está compuesto por una serie completa de observaciones equiespaciadas tomadas cada día en el período 1987-2000. Esto representa el formato requerido para aplicar DLNM en datos de series de tiempo.
3 Ejemplo 1: un DLM sencillo
En este primer ejemplo, especifico un DLM simple, evaluando el efecto de PM10 sobre la mortalidad, mientras se ajusta por el efecto de la temperatura. Para hacerlo, primero construyo dos matrices de base cruzada para los dos 2 predictores, y luego incluirlos en una fórmula modelo de una función de regresión. El efecto de PM10 se supone lineal en la dimensión del predictor, por lo que, desde este punto de vista, podemos definir esto como un DLM simple, incluso si el modelo de regresión también estima la función de retraso distribuida para la temperatura, que se incluye como un término no lineal.
Primero, crossbasis() para construir las dos matrices de
base cruzada, guardándolas en dos objetos. Los nombres de los dos
objetos deben ser diferentes para poder predecir las asociaciones por
separado para cada uno de ellos
cb1.pm <- crossbasis(chicagoNMMAPS$pm10, lag=15,
argvar=list(fun="lin"),
arglag=list(fun="poly",degree=4))cb1.temp <- crossbasis(chicagoNMMAPS$temp, lag=3,
argvar=list(df=5),
arglag=list(fun="strata",breaks=1))En aplicaciones con datos de series temporales, el primer argumento X se utiliza para especificar la serie de vectores.
La función pasa internamente los argumentos en argvary
arglag a onebasis()con el fin de construir la
base para predictor y retrasos, respectivamente. En este caso, asumimos
que el efecto de PM10 es lineal (fun="lin"), mientras
modela la relación con la temperatura a través de un spline cúbico
natural con 5 grados de libertad (fun="ns",elegido por
defecto). Los nudos internos (si no se proporcionan) son colocados por
ns() en los cuantiles equidistantes predeterminados,
mientras que los nudos límite se encuentran en el rango de temperatura,
por lo que solo df.debe epecificarse.
En cuanto a las bases para el espacio de los rezagos, especifico el
efecto rezagado de PM10hasta 15 días de retraso (mínimo retraso igual a
0 por defecto), con un función polinomial de grado 4 (configuración
degree=4).El efecto retardado de la temperatura está
definido por dos estratos de retardo (0 y 1-3), asumiendo que los
efectos son constantes dentro de cada estrato. El argumento
breaks = 1 define el límite inferior del segundo
intervalo.
La función de método proporciona una descripción general de las
especificaciones para la base cruzada (y las bases relacionadas en las
dos dimensiones).summary()para esta clase:
summary(cb1.pm)## CROSSBASIS FUNCTIONS
## observations: 5114
## range: -3.049835 to 356.1768
## lag period: 0 15
## total df: 5
##
## BASIS FOR VAR:
## fun: lin
## intercept: FALSE
##
## BASIS FOR LAG:
## fun: poly
## degree: 4
## scale: 15
## intercept: TRUEsummary(cb1.temp)## CROSSBASIS FUNCTIONS
## observations: 5114
## range: -26.66667 to 33.33333
## lag period: 0 3
## total df: 10
##
## BASIS FOR VAR:
## fun: ns
## knots: 0.2777778 6.666667 14.44444 20.94444
## intercept: FALSE
## Boundary.knots: -26.66667 33.33333
##
## BASIS FOR LAG:
## fun: strata
## df: 2
## breaks: 1
## ref: 1
## intercept: TRUEahora las dos bases cruzadas los objetos se pueden incluir en una
fórmula de modelo de un modelo de regresión. Los paquete
spline. En este caso, ajusté el modelo de serie temporal
asumiendo una distribución de Poisson sobredispersa, que incluye una
función de tiempo suave con 7 df/año (para corregir la estacionalidad y
la tendencia a largo plazo) y el día de la semana como factor:
library(splines)
model1 <- glm(death ~ cb1.pm + cb1.temp + ns(time, 7*14) + dow,
family=quasipoisson(), chicagoNMMAPS)La asociación estimada con niveles específicos de PM10sobre la
mortalidad, predicha por el modelo anterior, se puede resumir mediante
la crosspred()y guardado en un objeto con la misma
clase:
pred1.pm <- crosspred(cb1.pm, model1, at=0:20, bylag=0.2,
cumul=TRUE)La función incluye la basis1.pm y model1
objetos utilizados para estimar los parámetros como los dos primeros
argumentos, mientras que at= 0:20 establece que la
predicción debe calcularse para cada valor entero
de 0 a 20 $µgr/m^3$. Configurandobylag=0.2,la predicción se calcula a lo largo del espacio de retraso con un incremento de 0,2. Esta cuadrícula más fina está destinada a producir una curva de retraso más suave cuando se trazan los resultados. El argumentocumul(por defecto FALSO)indica que también se deben incluir asociaciones acumulativas incrementales a lo largo de los retrasos (nota: esta predicción solo se devuelve para retrasos enteros). No se define el centrado a través del argumento.cen,y el valor de referencia por lo tanto se establece en el valor 0 por defecto (esto sucede para la funciónlin()).Ahora que las predicciones se han almacenado enpred1.pm`,se
pueden trazar mediante funciones de método específicas. Por ejemplo:
plot(pred1.pm, "slices", var=10, col=3, ylab="RR",
ci.arg=list(density=15,lwd=2),
main="Lag-response curve for a 10-unit increase in PM10")
plot(pred1.pm, "slices", var=10, col=2, cumul=TRUE,
ylab="Cumulative RR",
main="Lag-response curve of incremental cumulative effects")
La función incluye pred1.pm objeto con los resultados
almacenados, y el argumento slices define que queremos
graficar la relación correspondiente a valores específicos de predictor
y retraso en el relacionado dimensiones.
Con var=10, muestra la relación de respuesta de retraso
para un valor específico de PM10, es decir 10gr/m3. Esta asociación se
define utilizando el valor de referencia de 0 gr/m3, proporcionando así
la asociación de predictor específico para un aumento de 10 unidades.
También elegí un color diferente para la primera trama. El argumento
cumul indica si las asociaciones acumulativas incrementale,
previamente guardadas en pred1.pm, debe ser ploteado. Los
resultados se muestran en las Figuras. Los intervalos de confianza se
establecen en el valor predeterminado valor “área” para el
argumento ci. En el panel izquierdo, los argumentos
adicionales se pasan al nivel bajo función de trazado
polygon() a través de ci.arg, para dibujar en
su lugar líneas de sombreado como intervalos de confianza.
La interpretación de estos gráficos es doble: la curva de retraso representa el aumento del riesgo en cada día futuro después de un aumento de PM_10 \(µgr/m^3\) en un día específico (interpretación directa), o en su defecto las aportaciones de cada día pasado con el mismo PM10 aumento del riesgo en un día específico (interpretación inversa).
Las parcelas en Figuras sugieren que el aumento inicial en el riesgo
de PM10 se invierte en retrasos más largos. El efecto acumulativo
general de un aumento de 10 unidades en PM10 más de 15 días de retraso
(es decir, sumandotodas las contribuciones hasta el retraso máximo),
junto con sus intervalos de confianza del 95% pueden ser extraídos
por los objetosallRRfit,allRRhighYallRRlowincluidos enpred1.pm,
`
pred1.pm$allRRfit["10"]## 10
## 0.9997563cbind(pred1.pm$allRRlow, pred1.pm$allRRhigh)["10",]## [1] 0.9916871 1.0078911example of application beyond time series - see vignette(“dlnmExtended”)
generate the matrix of exposure histories from the 5-year periods
Qnest <- t(apply(nested, 1, function(sub) exphist(rep(c(0,0,0,sub[5:14]), each=5), sub[“age”], lag=c(3,40))))
define the cross-basis
cbnest <- crossbasis(Qnest, lag=c(3,40), argvar=list(“bs”,degree=2,df=3), arglag=list(fun=“ns”,knots=c(10,30),intercept=FALSE)) summary(cbnest)
run the model and predict
library(survival) mnest <- clogit(case~cbnest+strata(riskset), nested) pnest <- crosspred(cbnest,mnest, cen=0, at=0:20*5)
bi-dimensional exposure-lag-response association
plot(pnest, zlab=“OR”, xlab=“Exposure”, ylab=“Lag (years)”) # lag-response curve for dose 60 plot(pnest, var=50, ylab=“OR for exposure 50”, xlab=“Lag (years)”, xlim=c(0,40)) # exposure-response curve for lag 10 plot(pnest, lag=5, ylab=“OR at lag 5”, xlab=“Exposure”, ylim=c(0.95,1.15))
example of extended predictions - see vignette(“dlnmExtended”)
compute exposure profiles and exposure history
expnested <- rep(c(10,0,13), c(5,5,10)) hist <- exphist(expnested, time=length(expnested), lag=c(3,40))
predict association with a specific exposure history
pnesthist <- crosspred(cbnest, mnest, cen=0, at=hist) with(pnesthist, c(allRRfit,allRRlow,allRRhigh))
example of user-defined functions - see vignette(“dlnmExtended”)
define a log function
mylog <- function(x) log(x+1)
define the cross-basis
cbnest2 <- crossbasis(Qnest, lag=c(3,40), argvar=list(“mylog”), arglag=list(fun=“ns”,knots=c(10,30),intercept=FALSE)) summary(cbnest2)
run the model and predict
mnest2 <- clogit(case~cbnest2+strata(riskset), nested) pnest2 <- crosspred(cbnest2, mnest2, cen=0, at=0:20*5)
plot and compare with previous fit
plot(pnest2, zlab=“OR”, xlab=“Exposure”, ylab=“Lag (years)”) plot(pnest2, var=50, ylab=“OR for exposure 50”, xlab=“Lag (years)”, xlim=c(0,40)) lines(pnest, var=50, lty=2) plot(pnest2, lag=5, ylab=“OR at lag 5”, xlab=“Exposure”, ylim=c(0.95,1.15)) lines(pnest, lag=5, lty=2)