I am using distributed lag non-linear models . I run a glm model with crossbasis matrix from DLNM package. When I am trying to get the predictions I get this error:
Error in crosspred(cbpm1, Tp1, by = 1, bylag = 1, at = speimin:speimax) : coef/vcov not consistent with basis matrix. See help(crosspred).
this happened when I tried for lag 1,2, and 3; but there were no error when i tried for lag 0, 4, and 5. I read the similar question from this link (https://stats.stackexchange.com/questions/482328/error-coef-vcov-not-consistent-with-basis-matrix). But still, I cannot figure it out with my own code. Your help is really meaningful for me. Thanks
The code is:
Dis <- ss$dis1
vkt <- equalknots(ss$T,nk=2)
lkt = logknots(1,nk=2)
vkpm <- equalknots(ss$spei3,nk=2)
lkpm <- logknots(1,nk=2)
speimin <- min(ss$spei3, na.rm = TRUE)
speimax <- max(ss$spei3, na.rm = TRUE)
cbt1 = crossbasis(ss$T, lag=1, argvar=list(fun="bs",degree=2,knots=vkt), arglag=list(knots=lkt))
cbpm1 <- crossbasis(ss$spei3, lag=1, argvar=list(fun="bs",degree=2,knots=vkpm), arglag=list(knots=lkpm))
Tp1 <- glm(Dis ~ cbt1 + cbpm1 + ns(RH,3)+ns(timeseries,2*5),
family=poisson(link=log),ss)
at=speimin:speimax
predsltp1 <- crosspred(cbpm1,Tp1,by=1,bylag=1,at=speimin:speimax)
here is the used library:
library(splines);library(class);library(stats);library(mda)
library(akima);library(gam);library(mgcv);library(foreign);library(som)
library(dlnm) #equalknots logknots crossbasis
library(splines) #ns
library(magrittr)
here is the reproducible sample of my dataset:
> a<-structure(list(job = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
> 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
> 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "all", class = "factor"),
> age3 = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
> 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
> 1L, 1L, 1L, 1L, 1L), .Label = "00_05", class = "factor"),
> sexA = structure(c(1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L,
> 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L,
> 2L, 1L, 2L, 2L, 1L), .Label = c("F", "M"), class = "factor"),
> All = c(65L, 53L, 92L, 68L, 81L, 103L, 144L, 92L, 44L, 40L,
> 54L, 19L, 55L, 61L, 72L, 89L, 77L, 68L, 71L, 27L, 15L, 18L,
> 39L, 52L, 52L, 58L, 27L, 44L, 32L, 37L), dis1 = c(6L, 0L,
> 9L, 0L, 0L, 0L, 9L, 0L, 3L, 6L, 3L, 0L, 0L, 3L, 6L, 0L, 9L,
> 3L, 0L, 3L, 6L, 0L, 0L, 0L, 0L, 0L, 3L, 0L, 0L, 0L), dis2 = c(3L,
> 6L, 0L, 0L, 0L, 0L, 0L, 3L, 0L, 0L, 0L, 3L, 0L, 0L, 6L, 6L,
> 0L, 0L, 0L, 3L, 3L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L),
> T = c(20.39032258, 20.39032258, 19.78387097, 19.78387097,
> 19.64193548, 19.64193548, 18.78709677, 18.78709677, 19.17419355,
> 19.17419355, 20.46774194, 21.63870968, 21.85806452, 21.85806452,
> 19.73448276, 19.73448276, 20.55357143, 20.55357143, 19.925,
> 29.12580645, 29.12580645, 29.39354839, 29.39354839, 28.96129032,
> 28.96129032, 27.36666667, 27.40333333, 27.40333333, 27.82333333,
> 27.82333333), RH = c(70.09677419, 70.09677419, 70.03225806,
> 70.03225806, 70.35483871, 70.35483871, 72.32258065, 72.32258065,
> 69.80645161, 69.80645161, 74.58064516, 77.58064516, 71.32258065,
> 71.32258065, 75.82758621, 75.82758621, 62.28571429, 62.28571429,
> 72.60714286, 77.61290323, 77.61290323, 75.06451613, 75.06451613,
> 75.61290323, 75.61290323, 76.03333333, 76.23333333, 76.23333333,
> 75.03333333, 75.03333333), PP = c(11.5, 11.5, 44.5, 44.5,
> 25.9, 25.9, 14, 14, 5, 5, 35.7, 34.1, 30.8, 30.8, 44.4, 44.4,
> 15.6, 15.6, 40.7, 184, 184, 137.1, 137.1, 377, 377, 110.5,
> 129.8, 129.8, 292, 292), spei3 = c(0.447495072, 0.447495072,
> 1.537295165, 1.537295165, 1.285067571, 1.285067571, 0.441010834,
> 0.441010834, 1.505630159, 1.505630159, 1.725831329, 1.075029338,
> -1.227673724, -1.227673724, 0.329690702, 0.329690702, 0.724314874,
> 0.724314874, 1.228544608, 0.60782059, 0.60782059, 0.191804009,
> 0.191804009, 1.752145476, 1.752145476, 1.94554333, 1.139058482,
> 1.139058482, -0.554472376, -0.554472376), timeseries = 1:30), class = "data.frame", row.names = c(NA,
> -30L))