MCMCglmm Poisson with an offset term and splines
Hello Jarrod,
Thank you so much for your patience!
This is the model that I used:
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)
I simply do not understand what am I doing wrong. These are my fixed terms:
1) age
2) x2
3) x8
4) x9
5) x3
6) l_lfvcspo
7) x4
8) x5
9) x6
10) x7
11) offset
There are thus 11 terms. If I add the intercept I obtain 12 fixed effects.
If I put k=12 per the number of fixed effects that I have, I get this error:
Error in MCMCglmm(y ~ age + x2 + x8 + x9 + x3 + l_lfvcspo + x4 + x5 + :
fixed effect mu prior is the wrong dimension
Do you think it might be an issue with the structure of the newdatab1 dataset? It seems Ok when I look at it. I do not know what else to do, I seem to have exhausted all options.
I apologize again, I really do not understand what am I doing wrong and how can I improve.
Best,
DNM
Sent from Outlook<http://aka.ms/weboutlook>
From: Jarrod Hadfield <j.hadfield at ed.ac.uk>
Sent: Wednesday, September 27, 2017 7:55 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
Sent: Wednesday, September 27, 2017 7:55 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, Yes - the intercept is a parameter so should be included. If you get 13 effects then the data/code you sent us are different from the data/code you are using. The issue is not really anything to do with your (mis)understanding about priors it is simply that you have set priors for k parameters whereas there are in fact not k parameters. Cheers, Jarrod On 27/09/2017 15:38, dani wrote: Hello again, Thanks! I was indeed referring to the k in the prior formula. I guess I am still confused whether I should add the intercept to the number of fixed effects, as I have 12 variables as I count as fixed effects without the intercept. Would it be possible to confirm that for me, please? Also, for k=12 in the prior formula, I get this error: Error in MCMCglmm(y ~ age + x2 + x8 + x9 + x3 + l_lfvcspo + x4 + x5 + : fixed effect mu prior is the wrong dimension I am puzzled about this error. I do not seem to understand how to select the fixed effects. Could someone point me to a source of information where I can learn more about this, please? Also, is there a book where I can learn more about prior specification? I think I need to see more examples. Thank you so much. I apologize for this, but I am very confused and I would like to learn more about this! Best, DNM Sent from Outlook<http://aka.ms/weboutlook> ________________________________ From: Jarrod Hadfield <j.hadfield at ed.ac.uk><mailto:j.hadfield at ed.ac.uk> Sent: Wednesday, September 27, 2017 7:28:17 AM To: dani; Ben Bolker Cc: Matthew; r-sig-mixed-models at r-project.org<mailto:r-sig-mixed-models at r-project.org> Subject: Re: [R-sig-ME] MCMCglmm Poisson with an offset term and splines 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><mailto: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<mailto: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><mailto: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<mailto: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"