[R-sig-ME] Mixed model with zero truncated Poisson distributed data

[Ben Bolker]

  [I'm taking the liberty of cc'ing this to the r-sig-mixed-models
mailing list, where it can be archived]

On 12-01-12 12:52 PM, Helen Ward wrote:
> Dear Ben,
> 
> Thank you for your speedy response. I've been having a go with MCMCglmm
> today, but I suspect I have over simplified the model (I'm new to
> Bayesian statistics), as I now get three suspiciously significant
> results...
> 
> The model I've tried (included below with its summary) runs, but my
> reading of the help file and the vignette make me suspect that I need to
> give some sort of prior of something to make it sensible. Am I on the
> right track?
> 
> I have also had a quick go with glmmADMB, trying the following models,
> 
>>  model1<-glmmadmb(RSperYr~Age+I(Age^2),data,family="truncpoiss",~1|DadID)
> 
> Error in switch(link, log = log, logit = qlogis, probit = qnorm, inverse
> = function(x) { :
> 
> EXPR must be a length 1 vector

   Here the problem is that you are specifying the random effect without
giving its name (i.e. random=~1|DadID), so R is trying to interpret it
as your desired link function (which is the fourth argument to the
glmmadmb() function: see ?glmmadmb).

> 
>>  model1<-glmmadmb(RSperYr~Age+I(Age^2)+(1|DadID),data,family="truncpoiss")
>>
> 
> Error in UseMethod("droplevels") :
> 
> no applicable method for 'droplevels' applied to an object of class
> "c('integer', 'numeric')",
> 
> but as you can see from the error mesddages I am obviously doing
> something a bit wrong here as well.

  As mentioned in another message on this list this morning, the problem
here is that you need to explicitly convert "DadID" to a factor from a
numeric value [e.g. data$DadID <- factor(data$DadID)] -- the current
error message is not very informative.

> I'll keep trying!, but are there any more obvious pointers you can give
> me please?
> 
> All the best,
> Helen
> 
> 
> This is the model I've tried using MCMCglmm.
>>
> model2<-MCMCglmm(RSperYr~Age+I(Age^2),random=~DadID,family="ztpoisson",data=data)
> 
>>  summary(model2)
> 
> Iterations = 3001:12991
> Thinning interval= 10
> Sample size= 1000
> 
> DIC: 808.5343
> 
> G-structure:~DadID
> 
> post.meanl-95% CI u-95% CI eff.samp
> 
> DadID0.07409 0.00057460.182239.75
> 
> R-structure:~units
> 
> post.mean l-95% CI u-95% CI eff.samp
> 
> units0.06205 0.0014050.176322.55
> 
> Location effects: RSperYr ~ Age + I(Age^2)
> 
> post.meanl-95% CIu-95% CI eff.samp pMCMC
> 
> (Intercept) -0.694625 -1.321090 -0.16384737.17 0.004 **
> 
> Age0.1978390.0667140.33041768.64 0.002 **
> 
> I(Age^2)-0.008292 -0.014887 -0.00158971.46 0.010 **
> 
> 

  At a *quick* glance this looks reasonable.  Try plot(model2$Sol) and
plot(model2$VCV) to see if the trace and density plots look reasonable
(compare them with the results from example("MCMCglmm"). You should take
a quick look at the Overview and CourseNotes vignettes (e.g.
vignette("Overview",package="MCMCglmm")) (the latter in particular is a
bit overwhelming but well worth digging into if you're going to use
these methods)

> 
> 
> 
> On 11/01/2012 18:46, Ben Bolker wrote:
>> Helen Ward<h.l.ward at ...>  writes:
>>
>>> I would like to describe the relationship between age and male
>>> reproductive success in a population of greater horseshoe bats.
>>>
>>> My data consists of three columns: MaleID, Age, NumberofPups (at that
>>> age). Many of the males appear multiple times in the data set, so I
>>> believe I need to derive a mixed model with MaleID as a random variable.
>>>
>>> The data is Poisson distributed, but zero-truncated. So far I have only
>>> succeeded in making a mixed model with a poisson distribution (using
>>> glmmPQL in the MASS package), and a zero truncated poisson model (using
>>> vglm in the VGAM package), but not a mixed model capable of handling
>>> zero truncated Poisson data.
>>>
>>> It has been suggested that I could just minus 1 from each value in the
>>> NumberofPups column to make a more usual Poisson distribution, so I can
>>> ignore the zero truncated bit. I have tried this and it changes the
>>> results of the model, but is this an acceptable transformation?
>>>
>>> If not, can anyone advise me on a mixed model that can handle zero
>>> truncated Poisson data please?
>>>
>>
>>    Thanks for letting us know about cross-posting.
>>
>>    You should be able to do this in either the MCMCglmm package or
>> (recent versions of) the glmmADMB package.  In MCMCglmm, use
>> family="ztpoisson"; in glmmADMB, use family="truncpoiss" ...
>>
>>    Ben Bolker
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>