[R-sig-ME] FW: Exponent random effect in nlmer
Ben Bolker
bbolker at gmail.com
Fri Nov 4 19:04:21 CET 2016
[cc'ing to r-sig-mixed-models]
OK, I'll put this back on the queue ... is there any chance that you
can send me/us a reproducible example?
The error message is driven from merPredD::updateDecomp in
src/predModule.cpp (line 291), in the package C++ code. This makes it
considerably harder to debug. There is no easily passable debug flag,
but if you download the source code and set debug=1 in line 258, you
will get at least a little bit more information about progress.
As for what the code is actually doing in the updateDecomp step
("update L, RZX, and RX"), your best bet is probably to look at
vignette("lmer",package="lme4"), eqs 48-52 and the code after it ... not
that it's easy ...
On 16-11-04 11:06 AM, Cole, Tim wrote:
> Dear Ben,
>
> I emailed you a while ago about this error message I was getting in
> nlmer: *Error in initializePtr() : Downdated VtV is not positive
> definite.* I'd be really grateful if you could respond - it's holding me
> up and I'd like to understand what is happening.
>
> Do you have any thoughts on what might be causing it? I can't find the
> code for initializePtr – is it accessible? Which matrix is VtV? Is it
> possible to set a breakpoint to track it?
>
> Any pointers would be really helpful…
>
> Best wishes,
> Tim
> ---
> Tim.cole at ucl.ac.uk <mailto:Tim.cole at ucl.ac.uk> Phone +44(0)20 7905 2666
> Fax +44(0)20 7905 2381
> Population, Policy and Practice Programme
> UCL Great Ormond Street Institute of Child Health, London WC1N 1EH, UK
>
> From: Tim Cole <tim.cole at ucl.ac.uk <mailto:tim.cole at ucl.ac.uk>>
> Date: Friday, 14 October 2016 14:41
> To: Ben Bolker <bbolker at gmail.com <mailto:bbolker at gmail.com>>
> Cc: "r-sig-mixed-models at r-project.org
> <mailto:r-sig-mixed-models at r-project.org>"
> <r-sig-mixed-models at r-project.org <mailto:r-sig-mixed-models at r-project.org>>
> Subject: Re: [R-sig-ME] Exponent random effect in nlmer
>
> Dear Ben,
>
> I have a model that can be written as:
>
> nlmer(y ~ (g - 1) ^ exp(k | id) , data=data, weights=w) (using a
> invalid but I hope obvious notation)
>
> where
> • -Inf < y < Inf
> • 0 < E(y) < 1
> • g is a grouping factor
> • k is a subject random effect
>
> k raises the group means to a subject-specific power while
> respecting the constraint on E(y).
>
> I've coded it in nlmer, but it gives *Error in initializePtr() :
> Downdated VtV is not positive definite*
>
> You have previously indicated that this can arise from badly
> rescaled parameters, so I'm wondering where I've gone wrong.
>
> Here is a simple example with 3 groups.
>
> > dim(data)
> [1] 759 4
> >
> > head(data)
> id y g w
> 1 1021 2.9968 01.02 0.2650
> 2 1022 0.8592 01.02 0.3931
> 3 1023 3.4657 01.02 0.2650
> 4 1024 0.3989 01.02 2.3648
> 5 1025 1.7230 01.02 0.2219
> 6 1026 0.7451 01.02 0.7046
> >
> > summary(moma <- with(data, model.matrix(~g - 1)))
> g01.02 g01.03 g02.03
> Min. :0.000 Min. :0.000 Min. :0.000
> 1st Qu.:0.000 1st Qu.:0.000 1st Qu.:0.000
> Median :0.000 Median :0.000 Median :0.000
> Mean :0.333 Mean :0.335 Mean :0.332
> 3rd Qu.:1.000 3rd Qu.:1.000 3rd Qu.:1.000
> Max. :1.000 Max. :1.000 Max. :1.000
> >
> > (start <- c(coef(lm(y ~g - 1, data=data, weights=w)), k=0))
> g01.02 g01.03 g02.03 k
> 0.7979 0.6355 0.9056 0.0000
> >
> > data <- cbind(data, moma)
> > nlmerfn <- function(g01.02,g01.03,g02.03, k) {
> + gc <- cbind(g01.02,g01.03,g02.03)
> + gc <- as.vector(as.matrix(gc * moma) %*% rep(1, ncol(moma)))
> + gc[gc <= 0] <- 1e-6
> + gk <- gc ^ exp(k)
> + grad <- moma * gk / gc * exp(k)
> + grad <- cbind(grad, k=gk * log(gk))
> + attr(gk, 'gradient') <- grad
> + gk
> + }
> > nlmer1 <- nlmer(y ~ nlmerfn(g01.02,g01.03,g02.03, k) ~
> + g01.02+g01.03+g02.03 + k + (k | id),
> data=data, weights=w, start=start)
> Error in initializePtr() : Downdated VtV is not positive definite
>
> Perhaps I need to include an intercept of some sort, but I'd very
> much value your thoughts.
>
> Best wishes,
> Tim
> ---
> Tim.cole at ucl.ac.uk <mailto:Tim.cole at ucl.ac.uk> Phone +44(0)20 7905
> 2666 Fax +44(0)20 7905 2381
> Population, Policy and Practice Programme
> UCL Great Ormond Street Institute of Child Health, London WC1N 1EH, UK
>
>
> From: Thierry Onkelinx <thierry.onkelinx at inbo.be
> <mailto:thierry.onkelinx at inbo.be>>
> Date: Wednesday, 12 October 2016 09:26
> To: Tim Cole <tim.cole at ucl.ac.uk <mailto:tim.cole at ucl.ac.uk>>
> Cc: "r-sig-mixed-models at r-project.org
> <mailto:r-sig-mixed-models at r-project.org>"
> <r-sig-mixed-models at r-project.org
> <mailto:r-sig-mixed-models at r-project.org>>
> Subject: Re: [R-sig-ME] Exponent random effect in nlmer
>
> Hi Tim,
>
> AFAIK nlmer requires the fixed and random effects to be
> additive. The model to be used _after_ this this summation can
> be non linear.
>
> Best regards,
>
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek / Research Institute for
> Nature and Forest
> team Biometrie & Kwaliteitszorg / team Biometrics & Quality
> Assurance
> Kliniekstraat 25
> 1070 Anderlecht
> Belgium
>
> To call in the statistician after the experiment is done may be
> no more than asking him to perform a post-mortem examination: he
> may be able to say what the experiment died of. ~ Sir Ronald
> Aylmer Fisher
> The plural of anecdote is not data. ~ Roger Brinner
> The combination of some data and an aching desire for an answer
> does not ensure that a reasonable answer can be extracted from a
> given body of data. ~ John Tukey
>
> 2016-10-11 12:30 GMT+02:00 Cole, Tim <tim.cole at ucl.ac.uk
> <mailto:tim.cole at ucl.ac.uk>>:
>
> Dear Thierry,
>
> Thanks very much for your speedy response.
>
> I agree my model looks odd, but it has a theoretical basis
> which I'd prefer not to spell out at this stage. Suffice to
> say that
> • -Inf < y < Inf
> • 0 < E(y) < 1
> • there is a subject random effect.
>
> For these reasons the usual models and/or transformations
> won't work, whereas my proposed exponent random effect ought
> to. I just need to fit it, to see if I'm right!
>
> Best wishes,
> Tim
> ---
> Tim.cole at ucl.ac.uk <mailto:Tim.cole at ucl.ac.uk> Phone
> +44(0)20 7905 2666 <tel:%2B44%280%2920%207905%202666> Fax
> +44(0)20 7905 2381 <tel:%2B44%280%2920%207905%202381>
> Population, Policy and Practice Programme
> UCL Great Ormond Street Institute of Child Health, London
> WC1N 1EH, UK
>
>
> From: Thierry Onkelinx <thierry.onkelinx at inbo.be
> <mailto:thierry.onkelinx at inbo.be>>
> Date: Tuesday, 11 October 2016 11:06
> To: Tim Cole <tim.cole at ucl.ac.uk <mailto:tim.cole at ucl.ac.uk>>
> Cc: "r-sig-mixed-models at r-project.org
> <mailto:r-sig-mixed-models at r-project.org>"
> <r-sig-mixed-models at r-project.org
> <mailto:r-sig-mixed-models at r-project.org>>
> Subject: Re: [R-sig-ME] Exponent random effect in nlmer
>
> Dear Tim,
>
> y centred on 0 and a valid range (0, 1) seems to be
> conflicting statements.
>
> Here a some solutions depending on y
>
> - y stems from a binomial process
> - use a binomial glmm.
> - y is continuous and you are willing to transform y
> - 0 < y < 1
> - apply a logit transformation on
> y. lmer(plogis(y) ~ f + (1 | id) )
> - 0 <= y < 1
> - apply a log transformation on y. lmer(log(y) ~
> f + (1 | id) )
> - 0 < y <= 1
> - apply a log transformation on 1 -
> y. lmer(log(1 - y) ~ f + (1 | id) )
> - y is continuous are not willing to transform y
> - use a beta regression with 0 and/or 1 inflation in
> case you have 0 or 1 in the data. Have a look at the
> gamlss package to fit this model.
>
> Best regards,
>
>
> ir. Thierry Onkelinx
> Instituut voor natuur- en bosonderzoek / Research
> Institute for Nature and Forest
> team Biometrie & Kwaliteitszorg / team Biometrics &
> Quality Assurance
> Kliniekstraat 25
> 1070 Anderlecht
> Belgium
>
> To call in the statistician after the experiment is done
> may be no more than asking him to perform a post-mortem
> examination: he may be able to say what the experiment
> died of. ~ Sir Ronald Aylmer Fisher
> The plural of anecdote is not data. ~ Roger Brinner
> The combination of some data and an aching desire for an
> answer does not ensure that a reasonable answer can be
> extracted from a given body of data. ~ John Tukey
>
> 2016-10-11 11:29 GMT+02:00 Cole, Tim <tim.cole at ucl.ac.uk
> <mailto:tim.cole at ucl.ac.uk>>:
>
> I have a model of the form
> m1 <- lmer(y ~ f + (1 | id) )
> where y is a continuous variable centred on zero, f
> is a unordered factor with coefficients b such 0 < b
> < 1, and there is a signficant random subject intercept.
>
> The random intercept can lead to predicted values
> outside the valid range (0, 1). For this reason I'd
> like to reformulate the model as
> m2 <- nlmer(y ~ (f - 1) ^ exp(1 | id) ) (using a
> invalid but I hope obvious notation), where the
> random effect is now a power centred on 1. This
> would constrain the fitted values to be within c(0, 1).
>
> My question is: can this be done in nlmer, and if so
> how? Please can someone point me in the right direction?
>
> Thanks,
> Tim Cole
> ---
> Tim.cole at ucl.ac.uk
> <mailto:Tim.cole at ucl.ac.uk><mailto:Tim.cole at ucl.ac.uk <mailto:Tim.cole at ucl.ac.uk>>
> Phone +44(0)20 7905 2666
> <tel:%2B44%280%2920%207905%202666> Fax +44(0)20 7905
> 2381 <tel:%2B44%280%2920%207905%202381>
> Population, Policy and Practice Programme
> UCL Great Ormond Street Institute of Child Health,
> London WC1N 1EH, UK
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org
> <mailto:R-sig-mixed-models at r-project.org> mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> <https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>
>
>
>
More information about the R-sig-mixed-models
mailing list