[R] Constrained non-linear mixed models
Spencer Graves
spencer.graves at pdf.com
Sun Jun 20 15:18:40 CEST 2004
WHY THIS FUNCTIONAL FORM?
If this were my problem, I would first want to know why we needed
this particular functional form and why (b*var3+a) and (b*var4+a) had to
be between 0 and 1. When I see things like this, my first impulse is
that these may only be approximations to something else, and then I
would seek a different approximation that might plausibly be better and
was not subject to the constraints.
CREATIVE PLOTS?
What kind(s) of random model(s) are you considering? What
kinds of plots have you made? Before I did heavy modeling, I'd compute
correlations and make plots designed to explore the chosen relationship
from different perspectives, e.g., using lattice graphics, contour and
perspective / wireframe plots, etc. Such plots might suggest
alternative functional forms.
PENALIZING "nlme"?
If I still needed to fit the model you suggested, I might first
rewrite it something like the following:
nlme(log(response)~z+y*log(var1)+x*log(var2)+log01(b*var3+a,
xlim)-x*log01(b*var4+a, xlim), ...
where log01 is something like the following [***UNTESTED***]:
log01 <- function(zz, xlim){
# log(zz) if xlim[1]<zz<xlim[2]
# else something that goes balistic othis range
# with rate controlled by xlim[3];
# log01 must be continuous and monotonic.
zz.sm <- (zz<xlim[1])
zz.lg <- (zz>xlim[2])
log.zz <- rep(NA, length(zz))
log.zz[zz.sm] <- (log(xlim[1])-
(exp(xlim[3]*(xlim[1]-zz[zz.sm]))-1))
log.zz[zz.lg] <- (log(xlim[2])+
(exp(xlim[3]*(zz[zz.lg]-xlim[2]))-1))
zz.gd <- !(zz.sm | zz.lg)
log.zz[zz.gd] <- log(zz[zz.gd])
log.zz
}
I would first run it with xlim something like c(0.1, 0.9, 0.01).
I'd then look at the cases for which I got log01 outside xlim[1:2] and
(gradually) adjust xlim to get something close to what I wanted, e.g.,
like c(0.01, 0.99, 100). If this didn't work well, I might also try to
modify log01 so it was also differentiable in zz.
FITTING SUBSETS?
I might also try fitting the same model to different subsets of
the data, perhaps suggested by the most plausible random effects models,
then use the parameter estimates as data for a second stage analysis,
following Box and Hunter, "A Useful Method for Model Building" paper;
you should be able to find a complete citation on the web.
MARKOV CHAIN MONTE CARLO (MCMC)?
I think the "gold standard" for this type of problem is MCMC.
There is software available for that, but I haven't used it and
therefore can't comment. However, before I got to that point, I would
make lots of plots and try to several alternative, simpler models to
convince myself that this was both appropriate and necessary.
hope this helps. spencer graves
Hwange Project wrote:
>Dear r-helpers,
>does anyone knows how to constrain the value of parameters in a non-linear
>mixed model:
>I've got something like that:
>nlme(log(response)~z+y*log(var1)+x*log(var2)+log((b*var3+a)/(b*var4+a)^x)
>where 0<x<1, a<1 and (b*(var3 or var4: range of values similar)+a) between 0
>and 1.
>or even better if you know how to linearize it !!
>
>and another one :
>What means ?
>Error in chol((value + t(value))/2) : the leading minor of order 1 is not
>positive definite
>(my statical books do not go far enough in the detail.
>It must be about the Cholesky something if I remember well what I read
>long-time ago).
>
>I'm sorry to bother you with this kind of things, but I'm in a remote place
>(~18°45'S,26°45'E)
>(even if e-mails are working at the incredible speed of 4.8 kbits/sec...)
>without a lot of statistical or computer ressources. Thanks,
>simon
>
>---------------------------------------------------------
>Simon Chamaillé
>Hwange Project
>CIRAD-CNRS
>Hwange Main Camp Research
>P. Bag 5776, Dete
>(+00 263) 18-647
>ciradhwg at mweb.co.zw <mailto:ciradhwg at mweb.co.zw>
>
>______________________________________________
>R-help at stat.math.ethz.ch mailing list
>https://www.stat.math.ethz.ch/mailman/listinfo/r-help
>PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
>
More information about the R-help
mailing list