convergence problem gamm / lme
Geert, Sorry for slow reply... I don't see any obvious problems with what you've done, so I guess it's the usual problem that PQL just doesn't *have* to converge, and the bit of extra flexibility of using a smooth is too much for it in this case. If you send me the data offline I can dig a little bit more if you like (I'll only use the data for this purpose etc. etc.) You are right that PQL does the same thing for Poisson and quasi-poisson. I don't think there is an easy way to use the values for the reduced dataset fit in the full dataset fitting, unfortunately. Another option is to use `gam' to fit the random effects. It'll be a bit slow with 70+ random effects, as you have, and it's a bit more work to set up, but it should converge. See ?gam.models which has some examples showing how to do this. best, Simon
On Thursday 29 January 2009 08:20, geert aarts wrote:
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=l ist(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=l ist(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.wood at bath.ac.uk To: r-help at 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(cod e_ 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 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 Sciences, University of Bath, Bath, BA2 7AY UK > +44 1225 386603 www.maths.bath.ac.uk/~sw283