[R-sig-ME] Exponent random effect in nlmer
Cole, Tim
tim.cole at ucl.ac.uk
Mon Nov 7 12:55:48 CET 2016
Thanks very much Ben. You've rather confirmed what I feared, that it's hairy to investigate.
I've gone back to my model and simplified it to be a weighted mean of 3 groups, with 10 subjects and either an additive or a multiplicative subject random effect. The additive model works fine while the multiplicative model fails with:
Error in initializePtr() : Downdated VtV is not positive definite
Have I got the multiplicative model right?
Best wishes,
Tim
> data
id y g w
1 1 1.71121 1 0.451575
2 2 1.02902 1 4.163295
3 3 0.80209 1 2.577572
4 4 0.62700 1 6.181288
5 5 0.50321 1 0.794809
6 6 1.08428 1 0.925991
7 7 0.86268 1 0.109342
8 8 0.59772 1 1.116947
9 9 0.24041 1 0.030605
10 10 2.73586 1 0.023148
11 1 2.38925 2 0.451575
12 2 0.89707 2 4.163295
13 3 0.56480 2 2.577572
14 4 0.82516 2 6.181288
15 5 -0.07298 2 0.794809
16 6 0.92487 2 0.925991
17 7 0.05092 2 0.109342
18 8 0.58145 2 1.116947
19 9 -1.50610 2 0.030605
20 10 5.32357 2 0.023148
21 1 1.39623 3 1.322324
22 2 0.87177 3 4.408423
23 3 0.70416 3 1.658296
24 4 1.31604 3 2.430068
25 5 -0.14503 3 0.201259
26 6 0.85298 3 1.088659
27 7 0.05903 3 0.081374
28 8 0.97278 3 0.399049
29 9 -6.26479 3 0.001769
30 10 1.94585 3 0.173258
>
> # additive subject random effect works
> lmer1 <- lmer(y ~g - 1 + (1|id), data=data, weights=w)
>
> # multiplicative subject random effect fails
> # model is effectively: y ~(g - 1) * (k|id) [where * indicates times not factor crossing]
>
> (start <- c(coef(lm1 <- lm(y ~g - 1, data=data, weights=w)), k=1))
g1 g2 g3 k
0.8084 0.7878 0.9923 1.0000
> moma <- with(data, model.matrix(~g - 1))
> data <- cbind(data, moma)
>
> nlmerfn <- function(g1,g2,g3,k) {
+ gk <- apply(cbind(g1,g2,g3) * moma, 1, sum) * k
+ attr(gk, 'gradient') <- cbind(moma * k, k=gk)
+ print(gk) # for monitoring
+ gk
+ }
>
> nlmer1 <- nlmer(y ~ nlmerfn(g1,g2,g3,k) ~ g1+g2+g3+k + (k|id), data=data, weights=w, start=start)
[1] 0.8084 0.8084 0.8084 0.8084 0.8084 0.8084 0.8084 0.8084 0.8084 0.8084 0.7878 0.7878 0.7878 0.7878 0.7878 0.7878 0.7878 0.7878 0.7878 0.7878 0.9923 0.9923 0.9923
[24] 0.9923 0.9923 0.9923 0.9923 0.9923 0.9923 0.9923
attr(,"gradient")
g1 g2 g3 k
1 1 0 0 0.8084
2 1 0 0 0.8084
3 1 0 0 0.8084
4 1 0 0 0.8084
5 1 0 0 0.8084
6 1 0 0 0.8084
7 1 0 0 0.8084
8 1 0 0 0.8084
9 1 0 0 0.8084
10 1 0 0 0.8084
11 0 1 0 0.7878
12 0 1 0 0.7878
13 0 1 0 0.7878
14 0 1 0 0.7878
15 0 1 0 0.7878
16 0 1 0 0.7878
17 0 1 0 0.7878
18 0 1 0 0.7878
19 0 1 0 0.7878
20 0 1 0 0.7878
21 0 0 1 0.9923
22 0 0 1 0.9923
23 0 0 1 0.9923
24 0 0 1 0.9923
25 0 0 1 0.9923
26 0 0 1 0.9923
27 0 0 1 0.9923
28 0 0 1 0.9923
29 0 0 1 0.9923
30 0 0 1 0.9923
Error in initializePtr() : Downdated VtV is not positive definite
>
> # multiplicative subject fixed effects for info
> data$ey <- fitted(lm1)
> (fe <- c(by(data, data$id, function(z) weighted.mean(z$y/z$ey, z$w))))
1 2 3 4 5 6 7 8 9 10
1.8810 1.0925 0.8193 0.9796 0.2187 1.1103 0.4286 0.7753 -0.9618 2.6167
---
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: Ben Bolker <bbolker at gmail.com<mailto:bbolker at gmail.com>>
Date: Friday, 4 November 2016 19:04
To: Tim Cole <tim.cole at ucl.ac.uk<mailto:tim.cole at ucl.ac.uk>>, "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: FW: [R-sig-ME] 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> <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> <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> <mailto:bbolker at gmail.com>>
Cc: "r-sig-mixed-models at r-project.org<mailto: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> <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> <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>
<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> <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>
<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>
<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>
<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> <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>
<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> <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>
<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>
<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>
<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 <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>
<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>
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list