[R-sig-ME] Exponent random effect in nlmer
Cole, Tim
tim.cole at ucl.ac.uk
Fri Oct 14 15:41:31 CEST 2016
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
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list