I am estimating a distributed lag non-linear model thanks to the R package dlnm. The independent variable is household consumption at time t, the independent variable is the weather conditions at time t-1 to t-L. The first step consists in using the matrix of weather exposure history to construct a cross-basis matrix where both the exposure and time dimensions can be modeled flexibly (crossbasis()). This cross-basis is plugged into a ols regression estimated with lm(). The dlnm package provides a few functions to plot and summarize the results (crosspred() and crossreduce()).
How do we get the parameters in the “original” scale, i.e. before the cross-basis transformation?
I found ways of doing it for a few specifications (see below), but is there a general method?
In the example below, the exposure dimension is modeled as a polynomial of degree 2 and the time dimension as a linear function.
library(dlnm)
cb <- crossbasis(chicagoNMMAPS$temp,lag=30,
argvar=list("poly",degree=2),
arglag=list("lin"))
model <- lm(cvd~cb,chicagoNMMAPS)
pred <- crosspred(cb,model,at=-20:30)
plot(pred,"slices",lag=0)
# my dirty way:
LAG<-0
SCALE<-attributes(cb)$argvar$scale
ce<-attributes(cb)$argvar$cen
B1<-(summary(model)$coeff[2,1]+summary(model)$coeff[4,1]*LAG)/SCALE
B2<-(summary(model)$coeff[3,1]+summary(model)$coeff[5,1]*LAG)/(SCALE^2)
B0<--(ce*B1+(ce^2)*B2)
xx<-(-20:30)
xx2<-xx^2
yhat<-B0+B1*xx+B2*xx2
lines(-20:30,yhat,lty=2,col="blue")
As shown on the graph below, the B's are correct (the blue-dotted line is perfectly aligned with the red line produced via the dlnm package). I found my “dirty way” after a lot of wrangling, not much mathematical sense behind it..
Is their a way to get the B's more directly and for all types of functional forms?
The solution is to set the scale parameter in the argvar function to 1. The crossreduce function then gives the parameters in the original scale. Below a an example:
# with a poly of degree 2 with natural splines at knots 3 and 10
library(dlnm)
cb <- crossbasis(chicagoNMMAPS$temp,lag=30,
argvar=list("poly",degree=2,scale=1),
arglag=list("ns",knots=c(3,10)))
model <- lm(cvd~cb,chicagoNMMAPS)
pred <- crosspred(cb,model,at=-20:30)
plot(pred,"slices",lag=3)
ce<-attributes(cb)$argvar$cen
xx<-(-20:30)
xx2<-xx^2
redvar<-crossreduce(cb, model,type="lag",value=3)
b1<-redvar$coefficients[1]
b2<-redvar$coefficients[2]
b0<--(ce*b1+(ce^2)*b2)
yhat<-b0+b1*xx+b2*xx2
lines(-20:30,yhat2,lty=2,col="green")