Simon, thanks for your reply and your suggestions.
I fitted the following glmm's
gamm3<-try(glmmPQL(count~offset(offsetter)+poly(lon,3)*poly(lat,3),random=list(code_tripnr=~1),family="poisson"))
Which worked OK (see summary below)
I also fitted a model using quasipoisson, but that didn't help. I actually also
thought that glmmPQL and gamm estimate the dispersion parameter and hence
assumes a quasipoisson distribution, even if you specify poisson. Is that
correct?
Finally I tried fitting a model to less data, and sometimes gamm managed to
converge (see summary below).
So would it be possible to use the parameter estimates from the model fitted to
less data as starting values for the gamm fitted to the full data set?
Or do you have any other suggestions?
Thanks.
Cheers Geert
>
gamm3<-try(glmmPQL(count~offset(offsetter)+poly(lon,3)*poly(lat,3),random=list(code_tripnr=~1),f
amily="poisson"))
iteration
1
iteration
2
iteration
3
> detach(Disc_age)
>
summary(gamm3)
Linear
mixed-effects model fit by maximum likelihood
Data: NULL
AIC BIC logLik
NA NA
NA
Random
effects:
Formula: ~1 | code_tripnr
(Intercept) Residual
StdDev:
0.001391914 231.9744
Variance
function:
Structure: fixed weights
Formula: ~invwt
Fixed
effects: count ~ offset(offsetter) + poly(lon, 3) * poly(lat, 3)
Value
Std.Error DF t-value p-value
(Intercept) -1.582 11.96 2024 -0.13232174 0.8947
poly(lon,
3)1 -4.048 1397.33 2024 -0.00289673 0.9977
poly(lon,
3)2 -22.013 699.71 2024 -0.03145996 0.9749
poly(lon,
3)3 -8.538 593.87 2024 -0.01437683 0.9885
poly(lat,
3)1 -109.624 666.05 2024 -0.16458856 0.8693
poly(lat,
3)2 -104.179 381.37 2024 -0.27316977 0.7848
poly(lat,
3)3 -10.661 221.93 2024 -0.04803585 0.9617
poly(lon,
3)1:poly(lat, 3)1 4290.737 61369.98 2024
0.06991589 0.9443
poly(lon,
3)2:poly(lat, 3)1 1853.559 36835.63 2024
0.05031972 0.9599
poly(lon,
3)3:poly(lat, 3)1 -240.521 25771.80 2024 -0.00933272 0.9926
poly(lon,
3)1:poly(lat, 3)2 2540.147 41378.38 2024
0.06138826 0.9511
poly(lon,
3)1:poly(lat, 3)2 2540.147 41378.38 2024
0.06138826 0.9511
poly(lon,
3)2:poly(lat, 3)2 -1803.911 21522.17
2024 -0.08381643 0.9332
poly(lon,
3)3:poly(lat, 3)2 1040.858 16352.56 2024
0.06365109 0.9493
poly(lon,
3)1:poly(lat, 3)3 632.587 12180.28 2024
0.05193535 0.9586
poly(lon,
3)2:poly(lat, 3)3 -394.339 13088.72 2024 -0.03012818 0.9760
poly(lon,
3)3:poly(lat, 3)3 -543.502 6221.71 2024 -0.08735569 0.9304
Correlation:
(Intr) ply(ln,3)1
ply(ln,3)2 ply(ln,3)3 ply(lt,3)1
poly(lon,
3)1 0.889
poly(lon,
3)2 0.938 0.878
poly(lon,
3)3 0.843 0.981
0.792
poly(lat,
3)1 -0.829 -0.949 -0.906
-0.882
poly(lat,
3)2 0.859 0.578 0.742
0.538 -0.474
poly(lat,
3)3 -0.552 -0.783 -0.579
-0.756 0.837
poly(lon,
3)1:poly(lat, 3)1 -0.947 -0.974
-0.940 -0.940 0.925
poly(lon,
3)2:poly(lat, 3)1 -0.934 -0.950
-0.857 -0.929 0.881
poly(lon,
3)3:poly(lat, 3)1 -0.818 -0.963
-0.866 -0.945 0.931
poly(lon,
3)1:poly(lat, 3)2 0.808 0.975
0.784 0.968 -0.928
poly(lon,
3)2:poly(lat, 3)2 0.737 0.575
0.853 0.465 -0.659
poly(lon,
3)3:poly(lat, 3)2 0.735 0.896
0.647 0.938 -0.765
poly(lon,
3)1:poly(lat, 3)3 -0.794 -0.592
-0.823 -0.518 0.591
poly(lon,
3)2:poly(lat, 3)3 -0.542 -0.737
-0.419 -0.781 0.635
poly(lon,
3)3:poly(lat, 3)3 -0.398 -0.383
-0.534 -0.334 0.425
ply(lt,3)2
ply(lt,3)3 p(,3)1:(,3)1 p(,3)2:(,3)1
poly(lon,
3)1
poly(lon,
3)2
poly(lon,
3)3
poly(lat,
3)1
poly(lat,
3)2
poly(lat,
3)3 -0.136
poly(lon,
3)1:poly(lat, 3)1 -0.708 0.690
poly(lon,
3)2:poly(lat, 3)1 -0.701 0.710 0.933
poly(lon,
3)3:poly(lat, 3)1 -0.499 0.738 0.956 0.849
poly(lon,
3)1:poly(lat, 3)2 0.458 -0.845
-0.915 -0.934
poly(lon,
3)2:poly(lat, 3)2 0.683 -0.344
-0.719 -0.522
poly(lon,
3)2:poly(lat, 3)2 0.683 -0.344
-0.719 -0.522
poly(lon,
3)3:poly(lat, 3)2 0.464 -0.655
-0.834 -0.884
poly(lon,
3)1:poly(lat, 3)3 -0.823 0.241 0.752 0.594
poly(lon,
3)2:poly(lat, 3)3 -0.300 0.707 0.612 0.788
poly(lon,
3)3:poly(lat, 3)3 -0.266 0.148 0.493 0.250
p(,3)3:(,3)1
p(,3)1:(,3)2 p(,3)2:(,3)2 p(,3)3:(,3)2
poly(lon,
3)1
poly(lon,
3)2
poly(lon,
3)3
poly(lat,
3)1
poly(lat,
3)2
poly(lat,
3)3
poly(lon,
3)1:poly(lat, 3)1
poly(lon,
3)2:poly(lat, 3)1
poly(lon,
3)3:poly(lat, 3)1
poly(lon,
3)1:poly(lat, 3)2 -0.928
poly(lon,
3)2:poly(lat, 3)2 -0.637 0.432
poly(lon,
3)3:poly(lat, 3)2 -0.851
0.935 0.245
poly(lon,
3)1:poly(lat, 3)3 0.642 -0.482 -0.894 -0.410
poly(lon,
3)2:poly(lat, 3)3 0.609 -0.822 0.007 -0.847
poly(lon,
3)3:poly(lat, 3)3 0.551 -0.327 -0.637 -0.291
p(,3)1:(,3)3
p(,3)2:(,3)3
poly(lon,
3)1
poly(lon,
3)2
poly(lon,
3)3
poly(lat,
3)1
poly(lat,
3)2
poly(lat,
3)3
poly(lon,
3)1:poly(lat, 3)1
poly(lon,
3)2:poly(lat, 3)1
poly(lon,
3)3:poly(lat, 3)1
poly(lon,
3)1:poly(lat, 3)2
poly(lon,
3)2:poly(lat, 3)2
poly(lon,
3)3:poly(lat, 3)2
poly(lon,
3)1:poly(lat, 3)3
poly(lon,
3)3:poly(lat, 3)1
poly(lon,
3)1:poly(lat, 3)2
poly(lon,
3)2:poly(lat, 3)2
poly(lon,
3)3:poly(lat, 3)2
poly(lon,
3)1:poly(lat, 3)3
poly(lon,
3)2:poly(lat, 3)3 0.080
poly(lon,
3)3:poly(lat, 3)3 0.684 -0.180
Standardized
Within-Group Residuals:
Min Q1 Med Q3 Max
-0.504980771 -0.000866948
0.028470924 0.078583094
33.247831244
Number
of Observations: 2113
Number
of Groups: 74
gamm3<-try(gamm(count~offset(offsetter)+s(lon,lat),random=list(code_tripnr=~1),family="quasipoisson",
niterPQL=200))
>
summary(gamm3$gam)
Family:
quasipoisson
Link
function: log
Formula:
count
~ offset(offsetter) + s(lon, lat)
Parametric
coefficients:
Estimate Std. Error t value Pr(>|t|)
X 1.31370
0.09854 13.33
>
summary(gamm3$lme)
Linear
mixed-effects model fit by maximum likelihood
Data: data
AIC
BIC logLik
2808.398 2837.845 -1398.199
Random
effects:
Formula: ~Xr.1 - 1 | g.1
Structure: pdIdnot
Xr.11 Xr.12
Xr.13 Xr.14 Xr.15
Xr.16 Xr.17 Xr.18
StdDev:
12.49623 12.49623 12.49623 12.49623 12.49623 12.49623 12.49623 12.49623
Xr.19 Xr.110
Xr.111 Xr.112 Xr.113
Xr.114 Xr.115 Xr.116
StdDev:
12.49623 12.49623 12.49623 12.49623 12.49623 12.49623 12.49623 12.49623
Xr.117 Xr.118
Xr.119 Xr.120 Xr.121
Xr.122 Xr.123 Xr.124
StdDev:
12.49623 12.49623 12.49623 12.49623 12.49623 12.49623 12.49623 12.49623
Xr.125 Xr.126
Xr.127
StdDev:
12.49623 12.49623 12.49623
Formula: ~1 | code_tripnr %in% g.1
(Intercept) Residual
StdDev: 0.8132693 5.077804
Variance
function:
Structure: fixed weights
Formula: ~invwt
Fixed
effects: list(fixed)
Value Std.Error
DF t-value p-value
XX 1.3137042 0.09863463 923
13.318894 0.0000
Xs(lon,lat)Fx1
-0.4406352 0.23114503 923 -1.906315
0.0569
Xs(lon,lat)Fx2
-0.6217519 0.24918031 923 -2.495189 0.0128
Correlation:
XX X(,)F1
Xs(lon,lat)Fx1 0.015
Xs(lon,lat)Fx2
-0.009 -0.148
Standardized
Within-Group Residuals:
Min Q1 Med Q3 Max
-3.42951750 -0.37448354
0.06432438 0.53690322 8.62026552
Number
of Observations: 1000
Number
of Groups:
g.1 code_tripnr %in% g.1
1 75
>
----------------------------------------
> From: s.w...@bath.ac.uk
> To: r-help@r-project.org
> Date: Fri, 23 Jan 2009 11:32:21 +0000
> Subject: Re: [R] convergence problem gamm / lme
>
> Geert,
>
> Can you get a simpler model with, say, a quadratic dependence on lon, lat to
> converge, using glmmPQL? The answer might give a clue about whether the issue
> is related to using a smoother, or is something more basic.
>
> How confident are you that the Poisson assumption is reasonable?
>
> Can the model be fitted to a random subsample of the data, or does it always
> fail? PQL can fail to converge, but it's usually not as obstinate as it seems
> to be in this case, if the model structure is reasonable and identifiable.
>
> best,
> Simon
>
>
>
>
>
> On Thursday 22 January 2009 15:52, geert aarts wrote:
>> Hope one of you could help with the following question/problem:
>> We would like to explain the spatial
>> distribution of juvenile fish. We have 2135 records, from 75 vessels
>> (code_tripnr) and 7 to 39 observations for each vessel, hence the random
>> effect for code_tripnr. The offset (�offsetter�) accounts for the haul
>> duration and sub sampling factor. There are no extreme outliers in lat/lon.
>> The model we try to fit is:
>>
>>
>>
>>
>> gamm3<-gamm(count~offset(offsetter)+s(lon,lat),random=list(code_tripnr=~1),
>>family="poisson", niterPQL=200)
>>
>> Maximum number of PQL iterations: 200
>>
>> iteration 1
>>
>> iteration 2
>>
>> Error in MEestimate(lmeSt, grps) :
>>
>>
>> NA/NaN/Inf in foreign function call (arg 1)
>>
>>
>>
>> We tried several things. We added some
>> noise to lon and lat, modelled the density instead of using a count with
>> model offset, and we normalized the explanatory variables. We also changed
>> several settings (see models below).
>>
>>
>>
>> Interestingly, we do manage to fit a more
>> complex model:
>>
>> gamm2<-gamm(count~offset(offsetter)+
>> s(lat,lon,year,dayofyear), random=list(code_tripnr=~1),family="poisson",
>> correlation = corGaus(0.1, form=~lat + lon))
>>
>>
>>
>> The models are fitted using mgcv 1.4-1 and
>> R 2.7.1 on a 64Bits Debian OS.
>>
>>
>>
>> So there seems to be a convergence problem, correct? And does someone have
>> an idea what might cause this? Secondly are there some tricks/solutions.
>> E.g. perhaps we could use the results from the more complex model (gamm2
>> above), but I do not know exactly how. All help/advice would be greatly
>> appreciated.
>>
>>
>>
>> Kind regards, Geert
>>
>>
>>
>>
>>
>>
>>
>> gamm3<-gamm(count~offset(offsetter)+s(lon,lat),
>> random=list(code_tripnr=~1),family="poisson", correlation = corExp(1,
>> form=~X + Y),nite
>>
>> rPQL=200)
>>
>> Maximum number of PQL iterations: 200
>>
>> iteration 1
>>
>> iteration 2
>>
>> Error in recalc.corSpatial(object[[i]],
>> conLin) :
>>
>>
>> NA/NaN/Inf in foreign function call (arg 1)
>>
>>> gamm3<-gamm(count~offset(offsetter)+s(lon,lat,k=c(1,1)),random=list(code_
>>>tripnr=~1),family="poisson",
>>
>> niterPQL=200)
>>
>> Maximum number of PQL iterations: 200
>>
>> iteration 1
>>
>> iteration 2
>>
>> Error in lme.formula(fixed = fixed, random
>> = random, data = data, correlation = correlation, :
>>
>> nlminb
>> problem, convergence error code = 1
>>
>>
>> message = false convergence (8)
>>
>> In addition: Warning messages:
>>
>> 1: In if (k < M + 1) { :
>>
>> the
>> condition has length> 1 and only the first element will be used
>>
>>
>>
>>
>>
>> .Options$mgcv.vc.logrange=0.001 # we also
>> tried higher settings
>>
>>
>> gamm3<-gamm(count~offset(offsetter)+s(lon,lat),random=list(code_tripnr=~1),
>>family="poisson", niterPQL=200, control=lmeControl(opt="optim"))
>>
>>
>>
>> Maximum number of PQL iterations: 200
>>
>> iteration 1
>>
>> iteration 2
>>
>> Error in optim(c(coef(lmeSt)),
>> function(lmePars) -logLik(lmeSt, lmePars),
>>
>>
>>
>> initial value in 'vmmin' is not finite
>>
>>
>>
>> gamm3<-gamm(count~offset(offsetter)+s(lon,lat),random=list(code_tripnr=~1),
>>family="poisson", niterPQL=200,control=lmeControl(minAbsParApV
>>
>> ar=0.0000000000001))
>>
>> Maximum number of PQL iterations: 200
>>
>> iteration 1
>>
>> iteration 2
>>
>> Error in recalc.corSpatial(object[[i]],
>> conLin) :
>>
>>
>> NA/NaN/Inf in foreign function call (arg 1)
>>
>>
>>
>>
>> gamm3<-gamm(count~offset(offsetter)+s(lon,lat),random=list(code_tripnr=~1),
>>family="poisson", niterPQL=200)
>>
>> Maximum number of PQL iterations: 200
>>
>> iteration 1
>>
>> iteration 2
>>
>> Error in MEestimate(lmeSt, grps) :
>>
>>
>> NA/NaN/Inf in foreign function call (arg 1)
>>
>>
>>
>>
>> gamm3<-gamm(count~offset(offsetter)+s(lon,lat,k=c(1,1)),random=list(code_tr
>>ipnr=~1),family="poisson", niterPQL=200)
>>
>> Maximum number of PQL iterations: 200
>>
>> iteration 1
>>
>> iteration 2
>>
>> Error in lme.formula(fixed = fixed, random
>> = random, data = data, correlation = correlation, :
>>
>>
>> nlminb problem, convergence
>> error code = 1
>>
>>
>> message = false convergence (8)
>>
>> In addition: Warning messages:
>>
>> 1: In if (k < M + 1) { :
>>
>> the
>> condition has length> 1 and only the first element will be used
>>
>> 2: In smooth.construct.tp.smooth.spec(object,
>> dk$data, dk$knots) :
>>
>>
>> basis dimension, k, increased to minimum possible
>>
>>
>>
>>
>>
>>
>> gamm3<-gamm(count~offset(offsetter)+s(lon,lat,k=c(8,8)),random=list(code_tr
>>ipnr=~1),family="poisson", niterPQL=200)
>>
>> Maximum number of PQL iterations: 200
>>
>> iteration 1
>>
>> iteration 2
>>
>> Error in lme.formula(fixed = fixed, random
>> = random, data = data, correlation = correlation, :
>>
>>
>> nlminb problem, convergence
>> error code = 1
>>
>>
>> message = false convergence (8)
>>
>> In addition: Warning messages:
>>
>> 1: In if (k < M + 1) { :
>>
>> the
>> condition has length> 1 and only the first element will be used
>>
>> 2: In 1:UZ.len : numerical expression has 2
>> elements: only the first used
>>
>> 3: In if (p.rank> ncol(XZ)) p.rank
>> <- ncol(XZ) :
>>
>> the
>> condition has length> 1 and only the first element will be used
>>
>> 4: In 1:p.rank : numerical expression has 2
>> elements: only the first used
>>
>> 5: In if (p.rank < k - j) Xf <- XZU[,
>> (p.rank + 1):(k - j), drop = FALSE] else Xf <- matrix(0, :
>>
>> the
>> condition has length> 1 and only the first element will be used
>>
>> 6: In (p.rank + 1):(k - j) :
>>
>>
>> numerical expression has 2 elements: only the first used
>>
>> 7: In 1:p.rank : numerical expression has 2
>> elements: only the first used
>>
>>
>>
>>
>> gamm3<-gamm(count~offset(offsetter)+s(lon,lat,k=c(4,4),fx=T),random=list(co
>>de_tripnr=~1),family="poisson", niterPQL=200)
>>
>> Maximum number of PQL iterations: 200
>>
>> iteration 1
>>
>> iteration 2
>>
>> Error in MEestimate(lmeSt, grps) :
>>
>>
>> NA/NaN/Inf in foreign function call (arg 1)
>>
>> In addition: Warning messages:
>>
>> 1: In if (k < M + 1) { :
>>
>> the
>> condition has length> 1 and only the first element will be used
>>
>> 2: In 1:UZ.len : numerical expression has 2
>> elements: only the first used
>>
>>
>>
>>
>> gamm3<-gamm(count~offset(offsetter)+te(lon,lat),random=list(code_tripnr=~1)
>>,family="poisson", niterPQL=200,control=lmeControl(opt="opti
>>
>> m"))
>>
>> Maximum number of PQL iterations: 200
>>
>> iteration 1
>>
>> iteration 2
>>
>> Error in optim(c(coef(lmeSt)),
>> function(lmePars) -logLik(lmeSt, lmePars),
>>
>>
>>
>> initial value in 'vmmin' is not finite
>>
>>
>>
>>
>> gamm3<-gamm(count~offset(offsetter)+s(lon,lat),random=list(code_tripnr=~1),
>>family="poisson", niterPQL=200,control=lmeControl(tolerance=
>>
>> 0.00000000000001))
>>
>> Maximum number of PQL iterations: 200
>>
>> iteration 1
>>
>> iteration 2
>>
>> Error in MEestimate(lmeSt, grps) :
>>
>>
>> NA/NaN/Inf in foreign function call (arg 1)
>>
>>> gamm3<-gamm(count~offset(offsetter)+s(lon,lat),random=list(code_tripnr=~1
>>>),family="poisson",
>>
>> niterPQL=200,control=lmeControl(niterEM=200))
>>
>> Maximum number of PQL iterations: 200
>>
>> iteration 1
>>
>> iteration 2
>>
>> Error in MEestimate(lmeSt, grps) :
>>
>>
>> NA/NaN/Inf in foreign function call (arg 1)
>>
>>
>>
>>
>> gamm3<-gamm(count~offset(offsetter)+s(lon,lat),random=list(code_tripnr=~1),
>>family="poisson",
>> niterPQL=200,control=lmeControl(msTol=0.00000000000000001))
>>
>> Maximum number of PQL iterations: 200
>>
>> iteration 1
>>
>> iteration 2
>>
>> Error in MEestimate(lmeSt, grps) :
>>
>>
>> NA/NaN/Inf in foreign function call (arg 1)
>>
>>
>>
>>
>> gamm3<-gamm(count~offset(offsetter)+s(lon,lat),random=list(code_tripnr=~1),
>>family="poisson",
>> niterPQL=200,control=lmeControl(.relStep=0.00000000000000000001))
>>
>> Maximum number of PQL iterations: 200
>>
>> iteration 1
>>
>> iteration 2
>>
>> Error in MEestimate(lmeSt, grps) :
>>
>>
>> NA/NaN/Inf in foreign function call (arg 1)
>>
>>
>>
>>
>> gamm3<-gamm(count~offset(offsetter)+s(lon,lat),random=list(code_tripnr=~1),
>>family="poisson",
>> niterPQL=200,control=lmeControl(nlmStepMax=0.00000000000000000001))
>>
>> Maximum number of PQL iterations: 200
>>
>> iteration 1
>>
>> iteration 2
>>
>> Error in MEestimate(lmeSt, grps) :
>>
>>
>> NA/NaN/Inf in foreign function call (arg 1)
>>
>>
>>
>>
>> gamm3<-gamm(count~offset(offsetter)+s(lon,lat),random=list(code_tripnr=~1),
>>family="poisson",
>> niterPQL=200,control=lmeControl(minAbsParApVar=0.0000000000001))
>>
>> Maximum number of PQL iterations: 200
>>
>> iteration 1
>>
>> iteration 2
>>
>> Error in MEestimate(lmeSt, grps) :
>>
>>
>> NA/NaN/Inf in foreign function call (arg 1)
>>
>>
>>
>>
>> gamm3<-gamm(count~offset(offsetter)+s(lon,lat),random=list(code_tripnr=~1),
>>family="poisson", niterPQL=200, control=lmeControl(returnObject=T))
>>
>> Maximum number of PQL iterations: 200
>>
>> iteration 1
>>
>> iteration 2
>>
>> Error in MEestimate(lmeSt, grps) :
>>
>>
>> Singularity in backsolve at level 0, block 1
>>
>> In addition: Warning messages:
>>
>> 1: In logLik.reStruct(object, conLin) :
>>
>>
>> Singular precision matrix in level -1, block 1
>>
>> 2: In logLik.reStruct(object, conLin) :
>>
>>
>> Singular precision matrix in level -1, block 1
>>
>> 3: In logLik.reStruct(object, conLin) :
>>
>>
>> Singular precision matrix in level -1, block 1
>>
>> 4: In logLik.reStruct(object, conLin) :
>>
>>
>> Singular precision matrix in level -1, block 1
>>
>> 5: In logLik.reStruct(object, conLin) :
>>
>>
>> Singular precision matrix in level -1, block 1
>>
>> 6: In MEestimate(lmeSt, grps) :
>>
>>
>> Singular precision matrix in level -1, block 1
>>
>>
>> _________________________________________________________________
>>
>>
>> [[alternative HTML version deleted]]
>
> --
>> Simon Wood, Mathematical Sciences, University of Bath, Bath, BA2 7AY UK
>> +44 1225 386603 www.maths.bath.ac.uk/~sw283
>
_________________________________________________________________
De leukste online filmpjes vind je op MSN Video!
http://video.msn.com/video.aspx?mkt=nl-nl
______________________________________________
R-help@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.
