Skip to content

MCMCglmm Poisson with an offset term and splines

9 messages · Jarrod Hadfield, dani, Matthew +1 more

#
Hello everyone,

I have a Poisson model with an offset term that involves repeated observations nested into two cross-classified groups.

The model includes
- four categorical variables
- 6 continuous variables (for one of them I would like to include a smoother)
- the offset=log(duration)

I first used the spl2 function to create the fixed ((x6numspline1) and random terms (x6numspline2) for the smoother. I added the random smoother term to the other two random intercepts (for student ID and classroom) that I have (which are cross-classified).

My question is: Do you find my model sound? Before I study the priors, I just wanted to run a default model - is my inclusion of an offset ok? Also, given that the observations are repeated and nested into both Student ID and classroom, I am not sure how to specify the variance structure in MCMCglmm (beginner here:))

mc_spl0 <- MCMCglmm(number_events ~ x1cat+x2cat+x8cat+x9cat+x3num+x4num+x5num+x6numspline1+x7num+x8num+log(duration),
                    random =~ ID+class+idv(x6numspline2),
                    data   = newdat,
                    family = "poisson",
                    thin   = 100,
                    burnin = 10000,
                    nitt   = 150000,
                    saveX=TRUE, saveZ=TRUE, saveXL=TRUE, pr=TRUE, pl=TRUE)

In addition, I am not sure what to make of the results for the offset term (included as a covariate in the model) in the output - how should I discuss them?

Thank you all so much!
Best regards,
N-M.
<http://aka.ms/weboutlook>
#
Hi,

The model looks OK as far as can be assessed without knowing the data. 
For the offset term you need to hold the associated coefficient at 1 by 
placing a strong prior on it. If  you want everything else to have the 
default prior then use:

k<-11 # number of fixed effects

prior<-list(B=list(V=diag(k)*1e8, 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$mu[k]<-1 # assuming the offset term is last
prior$B[k,k]<-1e-8

The interpretation of the offset is simply the coefficient is assumed to 
be one and that the rate at which events occur is constant.

Cheers,

Jarrod
On 22/09/2017 01:52, dani wrote:

  
    
#
Hi Jarrod,


Thank you so much for your prompt and helpful response!

I ran the code you sent me for the prior, but I am getting the following error:

Error in prior$B[k, k] <- 1e-08 :
  incorrect number of subscripts on matrix


Here is what I get for prior:
$B
$B$V
       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11]
 [1,] 1e+08 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
 [2,] 0e+00 1e+08 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
 [3,] 0e+00 0e+00 1e+08 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
 [4,] 0e+00 0e+00 0e+00 1e+08 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
 [5,] 0e+00 0e+00 0e+00 0e+00 1e+08 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
 [6,] 0e+00 0e+00 0e+00 0e+00 0e+00 1e+08 0e+00 0e+00 0e+00 0e+00 0e+00
 [7,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e+08 0e+00 0e+00 0e+00 0e+00
 [8,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e+08 0e+00 0e+00 0e+00
 [9,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e+08 0e+00 0e+00
[10,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e+08 0e+00
[11,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e+08

$B$mu
 [1] 0 0 0 0 0 0 0 0 0 0 0


$R
$R$V
[1] 1

$R$nu
[1] 0


$G
$G$G1
$G$G1$V
[1] 1

$G$G1$nu
[1] 0


$G$G2
$G$G2$V
[1] 1

$G$G2$nu
[1] 0


$G$G3
$G$G3$V
[1] 1

$G$G3$nu
[1] 0



$mu
 [1] NA NA NA NA NA NA NA NA NA NA  1



I am not sure what to do next. Any help would be very appreciated!

Best regards,

N-M



Sent from Outlook<http://aka.ms/weboutlook>
#
Hi Dani,

Jarrod probably meant

prior$B$mu[k]<-1 # assuming the offset term is last
prior$B$V[k,k]<-1e-8

Instead of:

prior$mu[k]<-1 # assuming the offset term is last
prior$B[k,k]<-1e-8


Sincerely,
Matthew

****************************************************
Matthew E. Wolak, Ph.D.
Assistant Professor
Department of Biological Sciences
Auburn University
306 Funchess Hall
Auburn, AL 36849, USA
Email: matthew.wolak at auburn.edu
On 22/09/17 12:01, dani wrote:
#
Hi Matt,

Thank you so much! I did try that code, as well, with the same result. Not sure what to do next:)

Best regards!


Sent from Outlook<http://aka.ms/weboutlook>
#
just a correction, the code prior$B$V[k,k]<-1e-8 works but the model provides this error


Error in MCMCglmm(y ~ x1 + x2 + X8 + x9 +x10 +  :
  prior list should contain elements R, G, and/or B only



Sent from Outlook<http://aka.ms/weboutlook>
#
Hi Dani,


Redo the prior specification from scratch (i.e., overwrite the entire 
`prior` object). The original code creates an unecessary/undefined 
element in the prior list.


Matthew

****************************************************
Matthew E. Wolak, Ph.D.
Assistant Professor
Department of Biological Sciences
Auburn University
306 Funchess Hall
Auburn, AL 36849, USA
Email: matthew.wolak at auburn.edu
On 22/09/17 12:17, dani wrote:

  
  
#
Hi again,


Thanks! It looks like it is working. This is like V looks like:
[,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] [,11]
 [1,] 1e+08 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
 [2,] 0e+00 1e+08 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
 [3,] 0e+00 0e+00 1e+08 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
 [4,] 0e+00 0e+00 0e+00 1e+08 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
 [5,] 0e+00 0e+00 0e+00 0e+00 1e+08 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00
 [6,] 0e+00 0e+00 0e+00 0e+00 0e+00 1e+08 0e+00 0e+00 0e+00 0e+00 0e+00
 [7,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e+08 0e+00 0e+00 0e+00 0e+00
 [8,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e+08 0e+00 0e+00 0e+00
 [9,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e+08 0e+00 0e+00
[10,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e+08 0e+00
[11,] 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 0e+00 1e-08

However, when I run my model, I get this new error:
 fixed effect V prior is not positive definite


Thanks!

Sent from Outlook<http://aka.ms/weboutlook>
#
I suspect that the broad range of the prior variances (from 1e8 to
1e-8) has led to an underflow error.  Try setting the prior variances
to 1e4 for all but the offset element and 1e-4 for the offset element ...

  Just a quick comment: this is why it can be helpful to provide a
minimal reproducible example
(you don't necessarily have to give us all of your data: see
https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example).
Debugging one step at a time can be frustrating for all concerned ...
On Fri, Sep 22, 2017 at 1:34 PM, dani <orchidn at live.com> wrote: