[R] Negative Binomial GAMM - theta values and convergence

[Simon Wood]

Hi Raeanne,

gamm fits using PQL, which doesn't always converge, however looking at 
your data, a couple of things stand out...

1. Trying a Tweedie seemed to give nicer residual plots than the 
negative binomial, and also converges with gamm. e.g.

SB.gam3<-gam(count~offset(vol_offset)+
   s(Depth_m, by=StnF, bs="cs")+StageF*RegionF+s(StnF,bs="re"),
              family=Tweedie(p=1.4),
              data=Zoop_2011[Zoop_2011$SpeciesF=='SB',],method="ML")
SB.gam3

plot(fitted(SB.gam3)^.5,residuals(SB.gam3))

qq.gam(SB.gam3)

... I just searched over p=1.2-1.6 in steps of .1 to find the optimum 
value for the Tweedie p parameter according to the ML score. see ?Tweedie

2. If you only need a random effect for StnF, and don't have to have the 
spatial autocorrelation structure then you can add the random effect in 
gam, as above - this will be much more efficient than using gamm 
(doesn't use PQL). See ?random.effects.

best,
Simon


On 07/12/12 12:00, Raeanne Miller wrote:
> Hi there,
>
> My question is about the 'theta' parameter in specification of a NB GAMM.
>
> I have fit a GAM with an optimum structure of:
>
> SB.gam4<-gam(count~offset(vol_offset)+
>    s(Depth_m, by=StnF, bs="cs")+StageF*RegionF,
>               family=negbin(1, link=log),
>               data=Zoop_2011[Zoop_2011$SpeciesF=='SB',])
>
> However, this GAM shows heterogeneity in the residuals, and the model fit is not great. I would like to add a random effect to this GAM, for station (StnF), and a spatial correlation structure, for Depth_m. When trying to fit an equivalent GAMM to start with, I used the formula:
>
>   SB.gamm1<-gamm(count~offset(vol_offset)+
>          s(Depth_m,by=StnF,bs="cs")+StageF*RegionF,
>                 family=negative.binomial(theta=1.41), niterPQL=50,
>                 data=Zoop_2011[Zoop_2011$SpeciesF=='SB',])
>
> The theta=1.41 came from SB.gam4, using the theta.ml() command. However, this model does not converge, even with 5000 iterations. When I increase theta to 6, the model converges, but I'm not sure why. I have also tried theta=c(1,10), which converges, but gives me lots of warnings:
>
> Warning messages:
> 1: In y + .Theta :
>    longer object length is not a multiple of shorter object length
>
>   Can anyone tell me why this might be?
>
> The model also does not ever converge if I replace family=negative.binomial()  with family=negbin(). I've looked at the help files for GAM and GAMM, and I'm not sure what the difference between these is?
>
> I have attached code to create my dataset, for those interested.
>
> With many thanks,
>
> Raeanne
> The Scottish Association for Marine Science (SAMS) is registered in Scotland as a Company Limited by Guarantee (SC009292) and is a registered charity (9206). SAMS has an actively trading wholly owned subsidiary company: SAMS Research Services Ltd a Limited Company (SC224404). All Companies in the group are registered in Scotland and share a registered office at Scottish Marine Institute, Oban Argyll PA37 1QA. The content of this message may contain personal views which are not the views of SAMS unless specifically stated. Please note that all email traffic is monitored for purposes of security and spam filtering. As such individual emails may be examined in more detail.
>
>
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>


-- 
Simon Wood, Mathematical Science, University of Bath BA2 7AY UK
+44 (0)1225 386603               http://people.bath.ac.uk/sw283