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.")}



-- 

*Andrew Crane-Droesch*
Energy and Resources Group
UC Berkeley
+1 215 435 2644
andrewcd@berkeley.edu
skype: andrew.crane-droesch
http://andrewcd.berkeley.edu



	[[alternative HTML version deleted]]

