R help - Dec 2012 - Negative Binomial GAMM - theta values and convergence

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.
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: Zoop_2011.txt
URL:
<https://stat.ethz.ch/pipermail/r-help/attachments/20121207/541f4f52/attachment-0002.txt>

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