[R-sig-ME] Singular estimated var-cov
Ken Beath
ken at kjbeath.com.au
Thu Oct 16 10:47:50 CEST 2008
On 16/10/2008, at 9:18 AM, francois.mercier at novartis.com wrote:
> 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.
>
I get slightly different results with the variance-covariance not
singular due to the time random effect variance being 0.0012559. This
is with lme4 version 0.999375-26.
The reason for the singularity problem is that the random effect
variance for time is close to zero, and doesn't need to be in the
model. If the model is fitted without it the log likelihood hardly
changes and the AIC and BIC are both lower.
Ken
> 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.
> <SerialData1.txt>_______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list