[R-sig-ME] Singular estimated var-cov
francois.mercier at novartis.com
francois.mercier at novartis.com
Thu Oct 16 00:18:07 CEST 2008
Thank you for your answers. I don't know why the attached file in my
initial post were not distributed. I attached hereafter the data set in
.txt format, hoping that it'll make it this time.
In addition, here are some further considerations on the problem I'm
facing with ...
1/ Load the data + run the original model which results in a 'singularity'
problem:
==================================================================
library(lme4)
dat <- read.csv("SerialData1.txt")
# Note: Random effects on both intercept and slope
lmer (formula = y ~ bas + drug + time + drug:time + (1+time|id), data =
dat)
Linear mixed-effects model fit by REML
Formula: y ~ bas + drug + time + drug:time + (1 + time | id)
Data: dat
AIC BIC logLik MLdeviance REMLdeviance
735.8 759.8 -359.9 714.7 719.8
Random effects:
Groups Name Variance Std.Dev. Corr
id (Intercept) 6.6569e+00 2.5801e+00
time 2.3614e-09 4.8594e-05 0.000
Residual 4.7228e+00 2.1732e+00
number of obs: 148, groups: id, 37
Fixed effects:
Estimate Std. Error t value
(Intercept) 3.55523 4.13153 0.861
bas 0.81279 0.15726 5.168
drug -1.53070 1.04169 -1.469
time -0.13033 0.10880 -1.198
drug:time 0.02914 0.15598 0.187
Correlation of Fixed Effects:
(Intr) bas drug time
bas -0.986
drug 0.163 -0.280
time -0.066 0.000 0.261
drug:time 0.046 0.000 -0.374 -0.697
Warning message:
In .local(x, ..., value) :
Estimated variance-covariance for factor ?id? is singular
2/ As per Gillian's suggestion, try out using centered time:
==================================================================
lmer (formula = y ~ bas + drug + ctime + drug:ctime + (1+ctime|id), data =
dat)
# ...estimated var-cov matrix is still singular
Linear mixed-effects model fit by REML
Formula: y ~ bas + drug + ctime + drug:ctime + (1 + ctime | id)
Data: dat
AIC BIC logLik MLdeviance REMLdeviance
735.8 759.8 -359.9 714.7 719.8
Random effects:
Groups Name Variance Std.Dev. Corr
id (Intercept) 6.6569e+00 2.5801e+00
ctime 2.3614e-09 4.8594e-05 0.000
Residual 4.7228e+00 2.1732e+00
number of obs: 148, groups: id, 37
Fixed effects:
Estimate Std. Error t value
(Intercept) 3.22942 4.12257 0.783
bas 0.81279 0.15726 5.168
drug -1.45786 0.96595 -1.509
ctime -0.13033 0.10880 -1.198
drug:ctime 0.02914 0.15598 0.187
Correlation of Fixed Effects:
(Intr) bas drug ctime
bas -0.988
drug 0.195 -0.302
ctime 0.000 0.000 0.000
drug:ctime 0.000 0.000 0.000 -0.697
Warning message:
In .local(x, ..., value) :
Estimated variance-covariance for factor ?id? is singular
3/ Is there a correlation between slopes and intercepts?
==================================================================
# fit linear models to individual patient data y=1+time
lmfit.1 <- lmList (y ~ 1 + time|id, data=dat)
# create vectors of intercepts and slopes
ints<-vector(length=37)
slps<-vector(length=37)
for (i in 1:37) ints[i]<-lmfit.1[[i]]$coefficients[1]
for (i in 1:37) slps[i]<-lmfit.1[[i]]$coefficients[2]
# scatterplot of slopes and intercepts by drug group
plot(ints,slps,col=2*(as.numeric(dat$drug[dat$time==0])+1))
lm0 <-
lm(slps[dat$drug[dat$time==0]==0]~ints[dat$drug[dat$time==0]==0])$coefficients
lm1 <-
lm(slps[dat$drug[dat$time==0]==1]~ints[dat$drug[dat$time==0]==1])$coefficients
curve(lm0[1]+lm0[2]*x,min(ints),max(ints),col=2,add=T)
curve(lm1[1]+lm1[2]*x,min(ints),max(ints),col=4,add=T)
# compute correlation coefficients
# cor for drug 0 is 0.2
cor(slps[dat$drug[dat$time==0]==0], ints[dat$drug[dat$time==0]==0])
# cor for drug 1 is -0.5
cor(slps[dat$drug[dat$time==0]==1], ints[dat$drug[dat$time==0]==1])
...is it possible that correlation ~0.5 as found below
is enough to produce a singular var-cov matrix?
4/ Random effect only on the slope
==================================================================
lmer (formula = y ~ bas + drug + time + drug:time + (time-1|id), data =
dat)
Linear mixed-effects model fit by REML
Formula: y ~ bas + drug + time + drug:time + (time - 1 | id)
Data: dat
AIC BIC logLik MLdeviance REMLdeviance
778.9 796.9 -383.5 761 766.9
Random effects:
Groups Name Variance Std.Dev.
id time 0.23548 0.48526
Residual 8.43043 2.90352
number of obs: 148, groups: id, 37
Fixed effects:
Estimate Std. Error t value
(Intercept) 3.93406 2.58167 1.524
bas 0.79816 0.09786 8.156
drug -1.50356 0.72968 -2.061
time -0.13033 0.18309 -0.712
drug:time 0.02914 0.26250 0.111
Correlation of Fixed Effects:
(Intr) bas drug time
bas -0.982
drug 0.115 -0.249
time -0.112 0.000 0.395
drug:time 0.078 0.000 -0.567 -0.697
# ...the fit is successful.
# note that correlation of fixed effects for "bas" and intercept is
almost 1
# does this mean that value of y before treatment ("bas") eliminates the
# need for a random effect on the intercept?
Best regards,
Francois
"Gillian Raab" <gillian.raab at googlemail.com>
10/13/2008 07:12 PM
To
"Douglas Bates" <bates at stat.wisc.edu>
cc
francois.mercier at novartis.com, r-sig-mixed-models at r-project.org
Subject
Re: [R-sig-ME] Singular estimated var-cov
I am sure the master (DB) is right in his comments But I wonder if you
have tried something easier. I did not get your data eitehr but you said
something about high correlations. Is your time variable centerd around
zero? This will make no difference to your fitted values and likelihood if
the models have converged to the ML solution, but centering will
sometimes overcome fitting problems.
If the variance covariance matrix at the solution really is singular then
obviously it won't help. If your fit gives you random effects than I would
suggest extracting them and plotting the equivalent lines, so you can see
what is going on. The most likley scenario from what you describe is that
your lines all pass through a common point for some choice of X value. You
could probably work this out from the fitted variance covariance matrix,
but I find the fitted lines help me to see what is going on better than
the algebra.
Gillian Raab, Edinburgh
On 13/10/2008, Douglas Bates <bates at stat.wisc.edu> wrote:
On Thu, Oct 9, 2008 at 7:27 AM, <francois.mercier at novartis.com> wrote:
> Dear list members,
> I try to fit a model (using lmer) to data recorded at 4 time points
> (days). Each such time series corresponds to a distinct subject. There
are
> two treatment groups. There is also a patient-level covariate ("o" or
> "b"). I am attaching the data frame (as a binary R object) and the R
> script that loads the data frame and fits the models.
I regret it has taken so long for you to get a response to your
question but I don't think that we can try the fit because you didn't
attach the data frame or the script - or at least they didn't make it
through the mail list software if you did include them.
> The questions are 1) whether the drug effect is influenced by the
> covariate, and 2) whether there is a temporal trend in drug effect over
> days.
> The problem is that according LMER the covariance matrix for this
problem
> is singular, and as a result the fitted models do not capture the
> variability of slopes that is seen in the data. Apparently there is a
> strong correlation between some parameters that leads to this
singularity
> ? Perhaps I misspecified the model for LMER (and LME) ?
It is possible for the estimated covariance matrix to be singular even
when there is significant variability in both the slope and the
intercept. An example of that is enclosed.
We can think of fitting mixed models as a smoothing problem where we
need to balance fidelity to the data against the complexity of the
model. The model complexity happens to be measured by a determinant
and a model with a singular covariance for the random effects has a
small value of this determinant. If there is not a correspondingly
large loss of fidelity to the data caused by the singular covariance
matrix then the estimates will be singular.
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
--
Gillian M Raab
10 Ainslie Place EH3 6AS
tel 0131 226 6234
mobile 07748 678 551
_________________________
CONFIDENTIALITY NOTICE
The information contained in this e-mail message is intended only for the
exclusive use of the individual or entity named above and may contain
information that is privileged, confidential or exempt from disclosure
under applicable law. If the reader of this message is not the intended
recipient, or the employee or agent responsible for delivery of the
message to the intended recipient, you are hereby notified that any
dissemination, distribution or copying of this communication is strictly
prohibited. If you have received this communication in error, please
notify the sender immediately by e-mail and delete the material from any
computer. Thank you.
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: SerialData1.txt
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20081015/a58f04b3/attachment.txt>
More information about the R-sig-mixed-models
mailing list