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}}
singular convergence in glmmPQL
2 messages · Austin, Peter, Bert Gunter
If I understand you correctly (your post was, ummm... a bit long), ?try should help. Condition on the class of the return inheriting from 'try-error' . i.e. result<-try(glimmPQL(...)) if(inherits(result,'try-error'))... As the try man page says, fancier things can be done via tryCatch() -- Bert Gunter Genentech Non-Clinical Statistics South San Francisco, CA "The business of the statistician is to catalyze the scientific learning process." - George E. P. Box
-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Austin, Peter
Sent: Monday, February 27, 2006 10:46 AM
To: r-help at stat.math.ethz.ch
Cc: Austin, Peter
Subject: [R] singular convergence in glmmPQL
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}}
______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html