[R-sig-ME] MCMCglmm Poisson with an offset term and splines
Jarrod Hadfield
j.hadfield at ed.ac.uk
Wed Sep 27 16:55:07 CEST 2017
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/8eff2170/attachment-0001.ksh>
More information about the R-sig-mixed-models
mailing list