Exponent random effect in nlmer
[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>