[R-sig-ME] MCMCglmm Poisson with an offset term and splines

[Jarrod Hadfield]

Hi,


There are 12 terms in the model not 13 so k=12. The k in the prior 
specification is completely unrelated to the number of knot points in 
the spline.


Cheers,


Jarrod


On 27/09/2017 15:16, dani wrote:
> Hello Jarrod,
>
> I have attached my code and my file.
>
> Variable age has 3 levels, variables x2, x8, and x9 have two levels 
> each, and the rest of the variables are continuous.
>
> I re-ran the spl2 function and this time around I obtained one single 
> fixed smoother (l_lfvcspo) and one single random smoother (l_lfvcspn). 
> Last time around I got 8 variables with suffixes from 1-8 for 
> l_lfvcspn and I did not know what to do with those.
>
> Also, as I had 11 fixed effects (corresponding to 11 variables), I 
> thought it was appropriate to choose k=11. I was not sure whether the 
> intercept needed to be counted as well as the levels of the 
> categorical fixed predictors (except for their reference categories). 
> Because I kept on getting the "mu" error for k=11, I tried k=13 (which 
> includes the intercept and the 2 levels for the 3-level age variable, 
> which I did not consider before). I am not sure this is the way to go, 
> I guess I need to read more to be able to model properly fixed effects 
> in R, but for now I was wondering whether this consideration of fixed 
> effects sounds ok in this particular example.
>
> The MCMC model worked properly based on the prior using k=13.
>
> Please see the code below:
>
> *#spl2 function*
>
> library(mgcv)
>
> spl2<-function(formula, data, p=TRUE, dataX=data){
> aug<-nrow(data)-nrow(dataX)
>   if(aug!=0){
>     if(aug<0){
>       stop("sorry nrow(dataX) must be less than or equal to nrow(data)")
>     }else{
>       augX<-matrix(0, aug, ncol(dataX))
> colnames(augX)<-colnames(dataX)
> dataX<-rbind(dataX, augX)
>     }
>   }
> smooth.spec.object<-interpret.gam(formula)$smooth.spec[[1]]
> sm<-smoothCon(smooth.spec.object, data=data, 
> knots=NULL,absorb.cons=TRUE, dataX=dataX)[[1]]
>   Sed<-eigen(sm$S[[1]])
>   Su<-Sed$vectors
>   Sd<-Sed$values
>   nonzeros <- which(Sd > sqrt(.Machine$double.eps))
>   if(p){
> Zn<-sm$X%*%Su[,nonzeros, drop=FALSE]%*%diag(1/sqrt(Sd[nonzeros]))
>   }else{
>     Zn<-sm$X[,-nonzeros, drop=FALSE]
>   }
> return(Zn[1:(nrow(data)-aug),,drop=FALSE])
> }
>
> *$spline terms*
> *
> *
> newdatab1$l_lfvcspo<-spl2(~s(f_lfv_c,k=10), data=newdatab1, p=F)
> newdatab1$l_lfvcspn<-spl2(~s(f_lfv_c,k=10), data=newdatab1)
> summary(newdatab1$l_lfvcspo)
> summary(newdatab1$l_lfvcspn)
>
> summary(newdatab1)
> dim(newdatab1)
> str(newdatab1)
>
> *#PRIOR*
>
>
> k<-13 # number of fixed effects
>
> prior<-list(B=list(V=diag(k)*1e4, mu=rep(0,k)),
>             R=list(V=1, nu=0),
> G=list(G1=list(V=1, nu=0),
>  G2=list(V=1, nu=0),
>  G3=list(V=1, nu=0)))
>
> prior$B$mu[k]<-1 # assuming the offset term is last
> prior$B$V[k,k]<-1e-4
>
> *# MCMC model*
>
> mcmc <- MCMCglmm(y ~ age+x2+x8+x9+x3+l_lfvcspo+x4+x5+x6+x7+offset,
>                  random =~ STUDYID+class+idv(l_lfvcspn),
>                  data   = newdatab1,
>                  family = "poisson", prior=prior)
> summary(mcmc)
>
>
> This model has worked, so I would like to thank you so much for all 
> your help so far. I guess I just want to make sure I understand how to 
> model the fixed effects in MCMCglmm and I also would like to make sure 
> that my model is correct.
>
> I truly appreciate all your invaluable help!
> Best regards,
> Dani NM
> ------------------------------------------------------------------------
> *From:* Jarrod Hadfield <j.hadfield at ed.ac.uk>
> *Sent:* Tuesday, September 26, 2017 12:01:19 PM
> *To:* dani; Ben Bolker
> *Cc:* Matthew; r-sig-mixed-models at r-project.org
> *Subject:* Re: [R-sig-ME] MCMCglmm Poisson with an offset term and 
> splines
>
> Hi Dani,
>
>
> It is still not possible for us to diagnose the problem. You need to 
> provide code+data that reproduces the error. f_lfv_c does not appear 
> in newdatab.
>
>
> Cheers,
>
>
> Jarrod
>
>
> On 25/09/2017 18:08, dani wrote:
>>
>> Hello again,
>>
>>
>> Thank you so much for your prompt response. I apologize for the silly 
>> questions, I am a true beginner and I am ashamed of my ignorance. I 
>> guess I should explain what I did:
>>
>>
>> I used the spl2 function and obtained the fixed and random factors 
>> corresponding to the variable I needed the smoother for (named f_lfv_c).
>>
>>
>> newdatab$l_lfvcspo<-spl2(~s(f_lfv_c,k=10), data=newdatab, p=F)
>> newdatab$l_lfvcspn<-spl2(~s(f_lfv_c,k=10), data=newdatab)
>>
>> I am not sure how to attach the random effects corresponding to the 
>> l_lfvcspn. I get this array of 8 variables and I am really not sure 
>> how to include them in the model. Should I get forget about spl2 and 
>> simply add the variable f_lfv_c as a fixed term and spl(f_lfv) in the 
>> idv random term?
>>
>> Also, it seems to me that I have 11 fixed effects, I am not sure what 
>> to do.
>>
>> I am really sorry about all these silly questions, I really do not 
>> understand how these things work, but I would like to know more about 
>> this.
>>
>> Best regards!
>>
>>
>> Sent from Outlook <http://aka.ms/weboutlook>
>> ------------------------------------------------------------------------
>> *From:* Jarrod Hadfield <j.hadfield at ed.ac.uk>
>> *Sent:* Monday, September 25, 2017 8:51:09 AM
>> *To:* dani; Ben Bolker
>> *Cc:* Matthew; r-sig-mixed-models at r-project.org
>> *Subject:* Re: [R-sig-ME] MCMCglmm Poisson with an offset term and 
>> splines
>>
>> Hi,
>>
>> The example is not reproducible: l_lfvcspn does not exist.
>>
>> The error is telling you that you don't have 11 fixed effects in the 
>> model. Change k to the number of fixed effects in the model.
>>
>> Jarrod
>>
>>
>>
>> On 25/09/2017 16:39, dani wrote:
>>> mc_spl1gna <- MCMCglmm(y ~ 
>>> age+x2+x8+x9+x3+l_lfvcspo+x4+x5+x6+x7+offset,
>>>        random =~ STUDYID+class+idv(l_lfvcspn),
>>>        data   = newdatab,
>>>        family = "poisson"
>>
>

-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20170927/803e3d48/attachment-0001.ksh>