Dear List,
I'm using gam in a multiple imputation framework -- specifying the knot
locations, and saving the results of multiple models, each of which is
fit with slightly different data (because some of it is predicted when
missing). In MI, coefficients from multiple models are averaged, as are
variance-covariance matrices. VCV's get an additional correction to
account for how variable they are between each other.
For this to work in the context of a penalized spline model, the knots
need to be specified identically for each model (this is assisted by
context knowledge), and each model needs to have the same number of
knots. This is what I've done, below. I run that code multiple times
with slightly different (imputed) datasets, but the number of
coefficients varies (between 263-265).
What gives? Why don't all of my models have the same number of
coefficients?
Thanks in advance!
Best,
Andrew
BCAR.knots = c(2,15,60,120)
INAR.knots = c(50,100,200,300)
bcph.knots = c(7.5,8.5,9.5,10.5)
htt.knots = c(350,450,550,650)
bc.prc.C.knots = c(.3,.45,.6,.8)
phi.knots = c(4.5,5.5,6.5,7.5)
CEC.knots = c(5,12,19,26)
soc.knots = c(10,20,30,40)
sand.knots = c(.2,.4,.6,.8)
clay.knots = c(.15,.3,.45,.6)
abslat.knots = c(10,20,30,45)
lon.knots = c(-50,0,50,125)
dum = as.vector(rep(1,length(trialid)))
doyee = NA
r.ints = tryCatch(gam(RR ~
as.factor(pot.trial)
+ as.factor(year)
+ as.factor(crop.legume)
+ as.factor(crop.fruit)
+ as.factor(feedstock)
+ s(trialid,bs="re",by=dum)
+ s(BCAR.imp,bs="cr",k=length(BCAR.knots),by=as.factor(pot.trial))
+
s(INAR.imp,bs="cr",k=length(INAR.knots),by=as.factor(crop.legume))
+ s(bcph.imp,bs="cr",k=length(bcph.knots),by=as.factor(year))
+ s(htt.imp,bs="cr",k=length(htt.knots))
+ s(bc.prc.C.imp,bs="cr",k=length(bc.prc.C.knots))
+ s(phi.imp,bs="cr",k=length(phi.knots),by=as.factor(year))
+ s(CEC.imp,bs="cr",k=length(CEC.knots))
+ s(soc.imp,bs="cr",k=length(soc.knots))
+ s(sand.imp,bs="cr",k=length(sand.knots),by=as.factor(year))
+ s(clay.imp,bs="cr",k=length(clay.knots))
+ s(abslat,bs="cr",k=length(abslat.knots))
+ te(INAR.imp,soc.imp,CEC.imp,k=c(4,4,4))
+ te(soc.imp,BCAR.imp,k=c(4,4))
+ te(phi.imp,bcph.imp,k=c(4,4))
+ te(clay.imp,CEC.imp,k=c(4,4))
+ te(htt.imp,bc.prc.C.imp,k=c(4,4),by=as.factor(feedstock))
,data=grain.data
,knots= list(
BCAR.imp=BCAR.knots,INAR.imp=INAR.knots,bcph.imp=bcph.knots,htt.imp=htt.knots,bc.prc.C.imp=bc.prc.C.knots,
phi.imp=phi.knots,CEC.imp=CEC.knots,soc.imp=soc.knots,sand.imp=sand.knots,clay.imp=clay.knots,
abslat.imp=abslat.knots
)
,method = "REML"
,weights = 1/nsame
), error=function(e) e, finally=doyee)
if(inherits(r.ints, "error")) {r.ints=doyee; print("an error happened
but it got handled.")}
