空间 GAM-GEE 中循环协变量的预测效果

如何解决空间 GAM-GEE 中循环协变量的预测效果

我在 {geepack} 包中安装了 GAM-GEE，因为我想在个体残差自相关中考虑一些已知个体的空间动物数据。我有一个循环协变量，我想使用循环样条拟合，因此我首先在 {mgcv} 中为样条生成基函数，因为 {splines} 包中目前不提供循环三次样条，然后将其用作 geeglm() 拟合中的协变量。我首先这样做了：

cyclic_basisfunc <- mgcv::gam(response ~s(cyclic_term,bs="cc"),fit=F,data=df,family = "binomial")$X %>%
  as_tibble() %>%
  dplyr::select(V2:ncol(.)) %>%
  as.matrix()

然后：

fit1 <- geepack::geeglm(response ~ cyclic_basisfunc +
                                   bs(beta1)+
                                   bs(lon+lat),family="binomial",id=as.factor(id),corstr = "independence")

我可以通过提取与该项对应的模型矩阵的相关元素（在本例中为循环项的 8 个基函数）来估计循环平滑的部分影响。但是，在使用 predict() 函数对空间网格进行预测时，我无法弄清楚如何适当地命名预测网格中的元素。循环项的输入变量是一个矩阵，任何将预测值包含为数字向量或矩阵的尝试都会返回错误。

我尝试在我的研究区域上生成一个预测网格（n.b. 尝试使用 cyclic_basisfunc 和 cyclic_term 的名称作为列标题）：

    predgrid <- data.frame(expand.grid(lat=seq(min(df$lat-0.1),max(df$lat+0.1),length.out = 20),lon=seq(min(df$lon-0.1),max(df$lon+0.1),cyclic_basis_func=seq(min(cyclic_basis_func),max(cyclic_basisfunc),length.out=20),beta1=25)))
predict(fit1,predgrid)

并得到错误：

Error: variable 'cyclic_basisfunc' was fitted with type "nmatrix.8" but type "numeric" was supplied

这是有道理的，因为原始协变量拟合为一个 8 列矩阵，每个基函数 1 个。已尝试将一列作为 8 列矩阵输入到预测网格中，但这似乎仍然不起作用。有没有办法为这个模型生成合适的预测网格？ （以下数据样本的 dput）。

df <- structure(list(response = c(0.117,0.915,0.22,0.284,0.551,0.871,0.25,0.261,0.117,0.875,0.67,0.533,0.324,0.138,0.154,0.286,0.935,0.744,0.118,0.224,0.865,0.13,0.889,0.115,0.703,0.917,0.812,0.14,0.219,0.114,0.24,0.163,0.228,0.106,0.243,0.908,0.472,0.95,0.217,0.133,0.92,0.26,0.958,0.113,0.147,0.132,0.496,0.148,0.119,0.186,0.721,0.162,0.157,0.269,0.129,0.357,0.188,0.361,0.137,0.128,0.872,0.121,0.135,0.466,0.15,0.589,0.134,0.122,0.463,0.369,0.813,0.145,0.911,0.17,0.123,0.649,0.476,0.367,0.923,0.895,0.923),cyclic_term = c(0.732222222222222,-2.28277777777778,2.56777777777778,-3.43333333333333,-0.89,-5.68611111111111,-3.47388888888889,-2.88277777777778,-1.79277777777778,1.50333333333333,-0.910555555555556,-4.14944444444444,-0.379027777777778,0.113333333333333,0.075,-3.99055555555556,-5.48388888888889,-1.84513888888889,0.286111111111111,-1.24833333333333,-2.19652777777778,-0.532222222222222,0.598333333333333,5.43222222222222,-2.73222222222222,-1.11125,3.16833333333333,-2.88055555555556,1.90319444444444,2.585,-5.64333333333333,-3.79666666666667,0.692083333333333,5.80666666666667,-4.42333333333333,1.95666666666667,2.61722222222222,-5.90055555555556,1.63,3.55222222222222,-4.20111111111111,-2.34,3.39277777777778,-4.32944444444444,1.32222222222222,4.74388888888889,0.251111111111111,0.0705555555555556,5.14305555555556,-3.92555555555556,-2.34277777777778,2.55777777777778,-3.75944444444444,2.32277777777778,1.82944444444444,-3.42111111111111,3.22388888888889,-3.90055555555556,1.94,-5.01888888888889,4.93902777777778,-2.97388888888889,4.11222222222222,1.55055555555556,-2.29055555555556,0.556666666666667,1.40375,3.52944444444444,4.56555555555556,1.30611111111111,-2.59944444444444,4.11166666666667,6.005,1.28111111111111,-2.35333333333333,-0.125,0.991666666666667,-4.29055555555556,4.64055555555556,1.19222222222222,-0.710555555555556,0.125,3.835,-3.79722222222222,1.46,0.455833333333333,-5.855,2.59277777777778,-1.42777777777778,4.815,0.508888888888889,-2.14333333333333,-1.47444444444444,5.01847222222222,3.06666666666667,5.92388888888889,1.39944444444444,5.00236111111111,4.21666666666667),beta1 = c(32.95,28.8,32.15,32.75,34.7,29.7,28.95,28.85,32.25,28.5,33.5,30.75,32.4,32.65,31.85,31.45,31.25,31.05,38.7,35.2,31.65,33,32.05,28.7,35.5,35.25,33.25,30.55,35.45,36,29.85,31.15,33.35,34.8,28.1,35.35,29.3,28.4,33.1,36.35,29,26.95,33.05,32.2,29.2,30.35,34.25,34.1,31.9,28.9,33.3,32.55,29.1,39.2,29.5,37.95,32.85,30.3,34.55,28.15,31.6,35.95,31.1,32.5,35.7,32.95,33.55),lon = c(-8.14386899769604,-8.1572279376935,-8.15157242384156,-8.11266145447895,-8.15174118001044,-8.15335952072546,-8.14600667978252,-8.15297882985764,-8.15568356665537,-8.13472008705139,-8.09368161533273,-8.10491923518749,-8.1603014305949,-8.15632063618518,-8.13543502792374,-8.09733172904193,-8.15868642071182,-8.12876868592058,-8.15690393084368,-8.10847025857788,-8.15564957894737,-8.16047965739412,-8.1538774272955,-8.13959002494812,-8.16031308174808,-8.13153629064039,-8.10225327088153,-8.11704735322503,-8.10320579591837,-8.09718723480212,-8.14769670963066,-8.15279598640478,-8.1536752924518,-8.15005347845117,-8.11959004402161,-8.15327362169542,-8.15338984397156,-8.15480377425017,-8.14843624352758,-8.1536150198704,-8.14265275265451,-8.15419676931013,-8.14388546959556,-8.12423783110794,-8.15865186565657,-8.12779523300791,-8.15498210353148,-8.15711005849511,-8.14498268712871,-8.12777905464012,-8.14658887045715,-8.14966988563538,-8.15137416124342,-8.14757286777317,-8.15659830243711,-8.11739216850327,-8.14816670318125,-8.15283383471808,-8.15503278645551,-8.12968355281506,-8.11532096238244,-8.15388445232111,-8.12550097443724,-8.14214966153336,-8.15406262640947,-8.15366204068896,-8.16073804747475,-8.14748077293754,-8.10236112317726,-8.13306401111526,-8.15754008293152,-8.11496173845736,-8.14744857712061,-8.13097980901942,-8.15712565747807,-8.16003438931358,-8.16002796870351,-8.11892837135011,-8.13700008392334,-8.15721941772399,-8.14819490909576,-8.1561399618782,-8.16012501716613,-8.11709369783742,-8.13470092070733,-8.14629675,-8.15713679162397,-8.13686372838595,-8.12430530737705,-8.15464372671311,-8.15669989585876,-8.16044796649062,-8.15766701538992,-8.09362511111111,-8.13870000839233,-8.14998125100998,-8.15243885077905,-8.12291705479452,-8.15384632745479),lat = c(33.6622974395803,33.6600173368609,33.6599819460086,33.6598656189251,33.6673908233704,33.6593042234088,33.6572338143965,33.6580708565473,33.6629485478539,33.6547317504883,33.6598712810572,33.6567040290043,33.6652851274788,33.6596383685524,33.6570077561196,33.6611995549193,33.661593106719,33.6588793953069,33.6614662478323,33.6584457075095,33.6641300278638,33.6621415089752,33.6598426484043,33.6570205688477,33.6727939605693,33.6593126847291,33.6591196082918,33.6591313969,33.6605346,33.6553558936485,33.6650662271625,33.6653484273653,33.6600826614748,33.6642524699626,33.6585006713867,33.658733988916,33.664975062454,33.6610514512704,33.6621421965042,33.6687249002733,33.6576841575676,33.6600360803426,33.6574947407709,33.6584060573279,33.6673802707071,33.6550894507841,33.6654066532252,33.6595139325288,33.6568044356436,33.6559013620991,33.6568076288124,33.6572189331055,33.6670448613725,33.6563744930416,33.6598780530983,33.6554991693199,33.6602992307323,33.6667773110577,33.6591215807971,33.6549378741848,33.6592876394984,33.6624833258972,33.6556429238423,33.6564236265265,33.6680011078708,33.6591718634038,33.6684195434343,33.6642363505652,33.6599223753418,33.6575357831156,33.6607284545898,33.6602992227631,33.6646009116676,33.6569601709497,33.6597380618501,33.6610439757783,33.6709440919706,33.6553216040898,33.6567993164063,33.6604019538554,33.6604042053223,33.6607238769535,33.6704235076904,33.6597095889516,33.6546409606935,33.658124,33.6659726122235,33.6550456763148,33.6591343852459,33.6648648385784,33.6633987426758,33.6695201510491,33.6609040792869,33.655895,33.6606597900391,33.6629373624763,33.6612642187822,33.6580512191781,33.6635307311956),id = c(8,14,12,7,8,10,6,4,2,9,5,3,11,12)),row.names = c("25101","15358","80097","89024","27479","98805","25425","86335","31333","93483","12171","100849","155418","12853","33858","100851","22470","149448","7443","12167","144934","33938","144132","91062","56909","153781","95039","99533","153760","32687","86913","12298","144133","7402","11672","13939","78548","98801","8135","24867","91818","32609","95688","11675","155218","94268","90367","7440","149447","145506","90571","35105","14210","26177","14975","16723","86359","34450","26139","14198","89237","145503","31062","92665","87694","78666","13917","155219","15350","60377","82820","33174","87056","7406","15370","15356","9330","33533","86726","95709","8131","96538","13911","14229","86539","93482","145837","22101","11305","144939","7391","7445","21817","32804","145314","98223","24175","8132","145504","87396"),class = "data.frame")

解决方法

您是否尝试过使用原始数据中使用的相同基础评估新值，然后将结果用作协变量？请注意，您会收到几个警告，因为您的响应是一个比例/概率，并且您没有给出权重或类似的东西。假设正确指定了模型，请查看下面的代码是否符合您的要求。

似乎您可以使用 predict.gam 中的“lpmatrix”类型从 GAM 中提取样条。看这里： How to extract fitted splines from a GAM (`mgcv::gam`) 老实说，我有点迷茫，所以我建议只为您的协变量设置基础而不拟合 GAM。

确保结的数量和位置有意义并建立基础。 据我了解，cSplinesDes() 在所需点评估基函数

spl_bs <- cSplineDes(x = df$cyclic_term,knots = seq(min(df$cyclic_term),max(df$cyclic_term),length.out = 8),ord = 4,derivs = 0)

取合理的名字

colnames(spl_bs) <- paste0("cyclic.",1:ncol(spl_bs))

您可能希望在数据框中包含这些新列。 请注意，我删除了下面的 Intercept 术语，但您可能希望将其保留在

df <- cbind(spl_bs,df)

将 GEE 与方便代码配合以指定型号

model <- reformulate(c(-1,paste0("cyclic.",1:(ncol(spl_bs))),"splines::bs(beta1)","splines::bs(lon+lat)"),response = "response")

fit1 <- geepack::geeglm(model,family="binomial",data=df,id=as.factor(id),corstr = "independence")

现在有了预测。创建新数据

predgrid <- data.frame(expand.grid(lat=seq(min(df$lat-0.1),max(df$lat+0.1),length.out = 20),lon=seq(min(df$lon-0.1),max(df$lon+0.1),cyclic_term=seq(min(df$cyclic_term),length.out=20),beta1=25))

设置与原始数据相同的基础，但现在评估新点

new_spline <- cSplineDes(predgrid$cyclic_term,derivs=0)
    
colnames(new_spline) <- paste0("cyclic.",1:ncol(new_spline))

组合新数据并根据 fit1 进行预测

newdf <- cbind(new_spline,predgrid)
predict(fit1,newdf)

注意纬度/经度超出了原始数据的范围，这是危险的