[R] singular convergence in glmmPQL

[Austin, Peter]

I am using the 'glmmPQL function in the 'MASS' library to fit a mixed effects logistic regression model to simulated data.  I am conducting a series of simulations, and with certain simulated datasets, estimation of the random effects logistic regression model unexpectedly terminates.  I receive the following error message from R:

Error in lme.formula(fixed=zz + arm.long,random=~1 | id.long, 
method="ML", : singular convergence (7)

I would like to modify my code so that the resultant model is assigned a null value, rather than having estimation terminate.

My R code for this example is contained below.

# Data are read in for the two sets of clusters.  There are 100 clusters
# in each of the two groups.

y.1j <-c(3,2,3,1,3,3,2,2,3,1,1,3,2,2,3,3,3,1,5,2,4,4,1,3,3,4,5,3,3,0,2,2,0,2,2,2,1,3,2,3,1,0,4,4,4,2,1,2,4,1,0,1,0,1,2,4,2,1,1,2,3,3,1,3,0,2,3,4,2,2,4,2,1,3,1,4,3,0,3,4,3,0,1,3,1,0,2,2,6,2,2,1,1,1,2,1,0,2,2,2)

y.2j <- c(3,4,3,3,3,4,1,8,2,3,3,1,3,2,1,2,4,6,1,3,6,3,4,3,5,5,4,6,3,3,4,0,4,2,2,4,1,0,5,2,7,4,4,3,3,4,4,6,1,1,2,2,2,4,0,2,1,5,5,2,3,4,1,3,1,1,1,4,3,3,3,2,1,3,1,3,2,3,2,3,4,2,3,2,7,2,2,2,3,2,6,2,2,3,3,3,2,3,1,4)
# Number of positive responses in each cluster in each group.

n.clusters.1 <- length(y.1j)
n.clusters.2 <- length(y.2j)
# Number of clusters in each group.

m.1j <- rep(20,n.clusters.1)
m.2j <- rep(20,n.clusters.2)
# 20 subjects in each cluster in each group.

M.1 <- sum(m.1j)
M.2 <- sum(m.2j)
# Number of subjects in each of the two groups.

K <- n.clusters.1 + n.clusters.2
# Total number of clusters in the two groups combined.

############################################################################
# Use a random effects logistic regression model with cluster-specific
# random intercepts. 
############################################################################

library("MASS")
# import MASS library for using function 'glmmPQL'

# create subject-specific responses from cluster-aggregate responses.
y.j <- c(y.1j,y.2j)
m.j <- c(m.1j,m.2j)
x.j <- m.j - y.j
resp.mat <- cbind(y.j,x.j)
# First column is number of successes, second column is number of failures.
# One row per cluster.

y <- rep(c(1,0),K)
pat.mat <- matrix(t(resp.mat),ncol=1,byrow=T)
y.long <- rep(y,pat.mat)
# Create vector of 1/0 outcomes with 1 observation per subject.

id.long <- rep(1:K,c(m.1j,m.2j))

arm.0 <- rep(1,M.1)
arm.1 <- rep(0,M.2)
arm.long <- c(arm.0,arm.1)
# Indicator variable for which group the cluster belongs to.

re.1 <- glmmPQL(y.long ~ arm.long,random = ~1|id.long,family=binomial)


My question is this:
Is it possible for the vector "re.1" to be returned with a null value if there is a singular convergence in the underlying call to 'lme'? Rather than having the execution terminate, can a null value be assigned to "re.1"?  As it stands, "re.1" is undefined when estimation terminates.

Thanks very much,

Peter

Peter Austin, PhD
Senior Scientist, Institute for Clinical Evaluative Sciences.
Associate Professor, Departments of Public Health Sciences and Health Policy, Management and Evaluation, University of Toronto.
 
Institute for Clinical Evaluative Sciences 
G106 - 2075 Bayview Avenue 
Toronto, Ontario, M4N 3M5 
Tel:  (416) 480-6131
Fax: (416) 480-6048 
ICES Web Site: www.ices.on.ca


___________________________________________________________________________ 
This email may contain confidential and/or privileged inform...{{dropped}}