5
$\begingroup$

I am performing distributed non-linear lag models in R.

I got the figure result of dlnm as shown in the vignette (pdf) on page 13:

enter image description here

The X-axis is lag, which I can understand. However, I cannot understand what the label for the Y-axis means. For the figure on page 13, the Y-axis is called "RR" (relative risk) which is positive; however, I got a negative RR when I input my data. My $X$ is temperature, my $Y$ is pollutant concentration. Can I still use the model?

Could anybody explain it to me a little bit?

The R code that I am using is:

OrigData <- read.csv("Hourly data for use SHULE-2.csv", na.strings = ".", 
                     header=T, sep=",")
CO2H6 = OrigData$CO2H6      # $
OutT  = OrigData$OutT       # $
Year  = OrigData$Year       # $
Month = OrigData$Month

argvar <- list(type="ns", df=11)
arglag <- list(type="ns", df=5)

suppressWarnings(cb1 <-crossbasis(OutT, lag=12, argvar=argvar, arglag=arglag))
model1 <- glm(CO2H6~ cb1 +as.factor(Month))
summary(model1)

cp <- crosspred(cb1,model1,-22.83:34.79,by=1, cumul=TRUE)
par(cex.axis=0.7,cex.lab=0.7)
d3 <- plot(cp, xlab="Delay Effects of OutT in IN2H-H6", zlab="CO2 con.", r=90, 
           d=0.3, col="red" , xlab="Out T")

plot(cp, "slices", var=10, ci="n", ylim= c(-1000,1000), ylab="CO2 con.", lwd=1.0, 
     xlim=c(0,10))
for(i in 4:6){
  lines(cp, "slices", var=c(-20,-10,0,10,20,30)[i], col=i+1, lwd=1.5)
}
legend(7,1000, col=5:7, lwd=0.9, pch=1:1, legend=c("OutT =10℃", "OutT =20℃", 
                                                   "OutT =30℃"), cex=0.7)

enter image description here

Since RR cannot be negative, could there be a way for me to transfer my x(Temperature) and y (Pollutant con.) to make it fitting the requirement of the model?

Thanks for all your help!!

$\endgroup$
  • 3
    $\begingroup$ Welcome. I went over this document: cran.r-project.org/web/packages/dlnm/vignettes/dlnmOverview.pdf, and on page 12, last paragraph it is written the explanation, but I did not understand too. Please, check if the 13's page figure (1a) is what you are talking about (if yes, you can use it as an example to your question). $\endgroup$ – Andre Silva Jun 25 '13 at 22:02
  • 5
    $\begingroup$ @user26221 your question is unfathomably cryptic. Andre has done you a service by at least identifying the package you're talking about. What function did you call? What figure do you mean? What does it look like? What are you attempting to achieve? $\endgroup$ – Glen_b -Reinstate Monica Jun 26 '13 at 2:25
  • 1
    $\begingroup$ @user26221 I've change the part concerning the interpretation of the graphics. I was unaware that the function crossbasis centeres the predictor at its mean and uses that as a reference value. To make it crystal clear: The y-axis is the predicted change of $\mathrm{CO}_{2}$ concentration for an increase of $z$ degrees compared to the mean temperature. $\endgroup$ – COOLSerdash Jun 26 '13 at 20:06
  • 1
    $\begingroup$ user26221 your question is looking improved, but there are still ways it doesn't seem to make sense. Your question talks about the relative risk being negative on your data, but I can't see any output that suggests that it was. What makes you say the RR is negative? $\endgroup$ – Glen_b -Reinstate Monica Jun 26 '13 at 21:16
  • 1
    $\begingroup$ If you do not set the error distribution family in glm the default is taken, which is a Gaussian distribution with the identity link function. This is equivalent with a normal linear regression (see my 4. point of my answer). So glm with no options is the same as lm. What family that you should take depends on the data. I think a Gaussian distribution is reasonable here. $\endgroup$ – COOLSerdash Jun 27 '13 at 17:23
8
$\begingroup$

Interpretation of the graph in your case

Note: The y-axis is not always the relative risk as in the example given in the vignette of the dlnm package. This is only the case in their example, because they used mortality data and Poisson regression models. In their framework, the exponentiated regression coefficient from the Poisson models $RR=\exp(\hat{\beta})$ is the relative risk. This is analogous to exponentiating the regression coefficients in logistic regression, which is the odds ratio.

Can I still use the model?

Yeah, you can still use such a model.

Let's summarize what you do:

  1. You use natural cubic B-splines as basis functions instead of polynomials to model the relationship between temperature and $\mathrm{CO}_{2}$ (arglag with option type="ns" instead of type="poly")
  2. You assume that the effect of temperature is non-linear, as you specify argvar as splines. One important thing you have to know for the interpretation of the plots is that the function crossbasis automatically centeres the values at the predictor mean (i.e. the mean temperature) if not specified otherwise. This is the reference value with which the predictions are later compared in the graphics.
  3. You consider lags up to 12 (option lag=12 in crossbasis). (Btw: Why do you suppress the warnings?)
  4. You calculate a GLM with Gaussian errors and the identity link function, which is equivalent to a simple linear regression (OLS). You could have used the lm function instead.

The plot that you have provided is interpreted as follows: The x-axis is the lag.

Interpretation of the values of the y-axis: The y-axis depicts the changes in $\mathrm{CO}_{2}$ concentration associated with an increase of 10, 20 or 30°C compared to the mean temperature. If predicted change is 0, this means that an increase in temperature is not associated with an increase in $\mathrm{CO}_{2}$ concentration compared to $\mathrm{CO}_{2}$ concentration at mean temperature: The predicted $\mathrm{CO}_{2}$ concentration is the same at $\bar{x}_{Temp}+z$ degrees (where $z$ is any amount, say 10 or 20 degrees) and at mean temperature $\bar{x}_{Temp}$.

This means that for an increase in temperature of 10°C, the temperature at lag 0 (on the same hour) increases the $\mathrm{CO}_{2}$ concentration compared to the mean temperature. Because you specified cumul=TRUE in crosspred, the effects are cumulative. The cumulative effects of an increase of 10°C are quasi nonexistent after 4 hours compared to the mean temperature. This suggests that the non-cumulative effects are negative at lags 1-4 and null effects from that on.

For temperature increases of 20 or 30°C, the cumulative effects on the $\mathrm{CO}_{2}$ concentration are lower in the first 1-4 hours compared to $\mathrm{CO}_{2}$ at the mean temperature. As with temperature increases of 10°C, the cumulative effects are practically nonexistent after 4 or 5 hours. Again: $\mathrm{CO}_{2}$ concentrations are the same at mean temperature and at an increase in temperature of 20 or 30°C after 4 or 5 hours.

I think a contour plot would be easier to interpret. Try the following code:

plot(cp, xlab="Temperature", col="red", zlab="CO2", shade=0.6,
     main="3D graph of temperature effect")

Interpretation of the example given in the vignette of the dlnm packge

First, a little something about distributed lag models. They have the form:

$$ Y_{i}=\alpha + \sum_{l=0}^{K}\beta_{j}x_{t-l} + \text{other predictors} +\epsilon_{i} $$ where $K$ is the maximum lag and $x$ is a predictor. This is just fitted using a multiple linear regression. So the coefficient $\beta_{1}$ would estimate the effect of $x_{t-1}$ of the day before on $Y_{t}$. In essence, multiple lags of the predictors are included in the model simultaneously. This obviously has the problem that the lagged predictors are highly correlated (autocorrelation).

A more advanced method are polynomial distributed lag models. It has the same basic formula as above, but the impulse-response function is forced to lie on a polynomial of degree $q$ (link to a paper for Stata):

$$ \beta_{i} = a_{0} + a_{1}i + a_{2}i^2 +\ldots+a_{q}i^q $$ where $q$ is the degree of the polynomial and $i$ the lag length. Another formulation is $$ \beta_{i} = a_{0} + \sum_{j=1}^{q}a_{j}f_{j}(i) $$

Where $f_{j}(i)$ is a polynomial of degree $j$ in the lag length $i$. A good introduction to the dlnm package and polynomial distributed lag models can be found here.

These models are often used in studies about air-pollution and health because air-pollution has lagged effects on health outcomes.

Let's look at this graph from the vignette of the dlnm package (page 13):

DLNM poylag

The degree of the polynomials was $q=4$ in this case so the green line is a polynomial of 4th degree. The y-axis is the relative risk (RR) estimated via Poisson regression and the x-axis the considered lag. The relative risk has the following interpretation: Persons who were exposed have a $(RR-1)\cdot100\%$ higher/lower chance of getting the outcome (e.g. death, lung cancer, etc.) compared to people who were not exposed. If $RR>1$ this means a positive association and if $RR<1$ means a protective association. A $RR=1$ means no association. We see that for every increase of $\textrm{PM}_{10}$ by 10 units ($\mu \mathrm{g}/m^{3}$), there is a $(1.001-1)\cdot100\%=0.1\%$ increase in the risk to die at lag 0 (i.e. on the same day as the exposure). Strangely, the exposure from about 9 days ago is protective: an increase of 10$~\mu \mathrm{g}/m^{3}$ is associated with a decreased risk to die compared to people with 10 units less exposure. We can also see that the exposure from 15 days before doesn't play a role (i.e. $RR\approx1$).

Let's look at the cumulative relative risk:

DLNM cumulative

This is the same as before but the effects are cumulated over time (i.e. summing all contributions from the lags up to the maximum lag). The red line starts at the same point as the green line in the first graphic (i.e. $\approx1.001$). We can see that people who have been exposed for five days have an increased cumulative risk of about $(1.005-1)\cdot100\%=0.5\%$ to die compared to non-exposed people. Because the green line goes below the relative risk of $1$ after a lag of about 5 days, the cumulative association after 15 days is nearly $1$. This means that the protective effects of $\textrm{PM}_{10}$ from lag 5 on have compensated the harmful effects from earlier lags. Whether that is scientifically reasonable is another quesiton.

$\endgroup$
  • 1
    $\begingroup$ +1 I hope that the OP can edit the question to match the high standard of the answer. $\endgroup$ – Glen_b -Reinstate Monica Jun 26 '13 at 7:55
  • $\begingroup$ @COOLSerdash Thank you so much for your exhausted explanation. I have read some paper in the field of air pollution and health and I saw the same figure as you uploaded: The RR is above 0 for their y is mortality, ranging from 0 to 1. For me, my y is concentration of pollutant (0-500ppm) and x is temperature. I got negative RR. That's really hard for me to interpretate it. eg, it is possitve at the first lag and goes down to negative at the second lag. The relationship between x and y should be negative based on reality. I have no idea how to interpretate it. $\endgroup$ – user26221 Jun 26 '13 at 15:03
  • $\begingroup$ @user26221 The relative risk can't be below 0. The y-axis is only the relative risk in this particular case, where Poisson regression models were done. In your case, it is certainly not the relative risk. To give you more advice, please edit your question and add your specific R code and the graphic as well. Then I can explain further. $\endgroup$ – COOLSerdash Jun 26 '13 at 15:25
  • $\begingroup$ @ COOLSerdash Thank you so much !!I will upload my R code immediately! $\endgroup$ – user26221 Jun 26 '13 at 15:33
  • $\begingroup$ @COOLSerdash Got it!I am sure I got a better understanding about it. But I don't know why the y shows negative value since all my input all possitive. I can't find an appropriate word to express my gratitude but I am really appreciate your exhausted explanation and kind help!! $\endgroup$ – user26221 Jun 26 '13 at 18:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.