[R] Constrained non-linear mixed models

[Spencer Graves]

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
>  
>