[R-sig-ME] MCMCglmm Poisson with an offset term and splines
Jarrod Hadfield
j.hadfield at ed.ac.uk
Wed Sep 27 17:19:13 CEST 2017
Hi,
If I use the code you sent me with the data you sent me there are 12
effects and using k=12 works fine. You must have changed the data and/or
the code.
Jarrod
On 27/09/2017 16:13, dani wrote:
>
>
> 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
>
> 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>
>> *Sent:* Wednesday, September 27, 2017 7:28:17 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,
>>
>>
>> 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/22d1115d/attachment-0001.ksh>
More information about the R-sig-mixed-models
mailing list