Skip to content
Prev 8446 / 29559 Next
Original Raw

Help with glmmPQL using Scottish Lip Cancer data

Corey Sparks 
Dear List,
I'm trying to replicate results from several sources that analyze the
Clayton and Kaldor (1987) Scottish Lip Cancer data using a Poisson GLMM with
spatial correlation structure using a simple exponential variogram model.
The only way (I think) that my data are different from the published data
are that i'm using the Lat/Long instead of the OSGB projection.

Here is what I've tried:
scot<-read.csv("scotland.csv")
names(scot)
scot$aff<-scot$AFF/10
scot$x<-scot$Longitude
scot$y<-scot$Latitude
scot$smr<-scot$Observed/scot$Expected
scot$logsmr<-log(scot$smr+1)
coordinates(scot)=~x+y

library(MASS)
library(nlme)

#This is the first try, using the corEXP() directly in glmmPQL

fit<-glmmPQL(Observed~offset(log(Expected))+aff,family=poisson, data=scot,
random=~1|District, correlation=corExp(form=~x+y), niter=200)

#it fails with:
iteration 1
Error in corFactor.corSpatial(object) :
  NA/NaN/Inf in foreign function call (arg 1)
In addition: Warning messages:
1: In min(unlist(attr(object, "covariate"))) :
  no non-missing arguments to min; returning Inf
2: In min(unlist(attr(object, "covariate"))) :
  no non-missing arguments to min; returning Inf


#Then I try to initialize the spatial correlation as in ASDAR
sp2<-corSpatial( form=~x+y, type="exponential")
spcor<-Initialize(sp2, as(scot, "data.frame")[,c("x","y")])

glmmPQL(Observed~offset(log(Expected))+aff,family=poisson, data=scot,
random=~1|District, correlation=spcor, niter=200, na.action="na.omit")

#This fails with:
iteration 1
iteration 2
iteration 3
iteration 4
iteration 5
Error in lme.formula(fixed = zz ~ offset(log(Expected)) + aff, random = ~1 |
: 
  nlminb problem, convergence error code = 1
  message = function evaluation limit reached without convergence (9)

Basically, i'm wondering if  anyone can get this to work, because i'm at my
wits end.  


My system info is:
sessionInfo()
R version 2.11.1 (2010-05-31)
x86_64-apple-darwin9.8.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] nlme_3.1-96  MASS_7.3-6   gstat_0.9-69 sp_0.9-64

loaded via a namespace (and not attached):
[1] grid_2.11.1    lattice_0.18-8 tools_2.11.1


Best, 
Corey

I'm pasting in the data, in case anyone wants to have a look:

District,Observed,Expected,AFF,Latitude,Longitude
1,9,1.4,16,57.29,5.5
2,39,8.7,16,57.56,2.36
3,11,3,10,58.44,3.9
4,9,2.5,24,55.76,2.4
5,15,4.3,10,57.71,5.09
6,8,2.4,24,59.13,3.25
7,26,8.1,10,57.47,3.3
8,7,2.3,7,60.24,1.43
9,6,2,7,56.9,5.42
10,20,6.6,16,57.24,2.6
11,13,4.4,7,58.12,6.8
12,5,1.8,16,58.06,4.64
13,3,1.1,10,57.47,3.98
14,8,3.3,24,54.94,5
15,17,7.8,7,56.3,3.1
16,9,4.6,16,57,3
17,2,1.1,10,57.06,4.09
18,7,4.2,7,55.65,2.88
19,9,5.5,7,57.24,4.73
20,7,4.4,10,55.35,2.9
21,16,10.5,7,56.75,2.98
22,31,22.7,16,57.12,2.2
23,11,8.8,10,56.4,5.27
24,7,5.6,7,55.63,3.96
25,19,15.5,1,56.2,3.3
26,15,12.5,1,56.1,3.6
27,7,6,7,55.24,4.09
28,10,9,7,55.95,2.8
29,16,14.4,10,56.6,4.09
30,11,10.2,10,55.9,3.8
31,5,4.8,7,55.47,4.55
32,3,2.9,24,55,4.36
33,7,7,10,55.83,3.2
34,8,8.5,7,56.3,4.73
35,11,12.3,7,55.29,4.98
36,9,10.1,0,55.94,4.95
37,11,12.7,10,55.76,5.02
38,8,9.4,1,55.91,4.18
39,6,7.2,16,56.15,4.99
40,4,5.3,0,56.05,4.91
41,10,18.8,1,55.88,4.82
42,8,15.8,16,56.03,4
43,2,4.3,16,56.15,3.96
44,6,14.6,0,55.82,4.09
45,19,50.7,1,55.93,3.4
46,3,8.2,7,55.65,4.75
47,2,5.6,1,55.71,4.45
48,3,9.3,1,55.79,4.27
49,28,88.7,0,55.9,4.55
50,6,19.6,1,56.45,3.2
51,1,3.4,1,56,4.27
52,1,3.6,0,56.15,4.64
53,1,5.7,1,55.79,4.7
54,1,7,1,55.99,4.45
55,0,4.2,16,55.68,3.38
56,0,1.8,10,55.18,3.4