[R-sig-ME] GLMM with bounded parameters

[Pierre de Villemereuil]

Dear all, 

Thank you for your interesting suggestions!

I finally resorted to JAGS to analyse the data, which gives me more control on 
the estimation process (e.g. distinguishing between location and shape 
parameters).

Here's the BUGS model for people interested:

data {
    N <- length(Y)
}

model{
    for (i in 1:N) {
        lat[i] <- scale*exp(-((X[i] - location)/sigma)^2)
        Y[i] ~ dpois(lat[i])
    }
    
    scale ~ dunif(0,20)
    sigma ~ dunif(0,2)
    location ~ dunif(-2,2)
}

Cheers,
Pierre

Le lundi 19 décembre 2016, 10:05:55 CET François Rousset a écrit :
> Le 16/12/2016 à 13:43, Pierre de Villemereuil a écrit :
> > Dear all,
> > 
> > I'm trying to fit a predictive bell curve on count data with Poisson noise
> > around the curve. The idea is to estimate the optimum of fitness according
> > to some trait.
> > 
> > The good news is that I don't need to resort to non-linear modelling to do
> > that, because the exponential link combined with a polynomial formula can
> > do the job as shown in the dummy example below:
> > x <- runif(300,0,10)
> > fitness <- rpois(300,10*dnorm(x,mean=5,sd=2))
> > 
> > mod = glm(fitness ~ x + I(x^2),family="poisson")
> > plot(fitness ~ x)
> > points(predict(mod,type="response") ~ x,col="red")
> > 
> > My problem is that I'd like to impose some constraints on the parameters
> > to
> > ensure that a bell-shape is fitted, e.g. that the parameter for x is
> > positive and the parameter for I(x^2) is negative.
> 
> One can use offsets and a brute-force optim() over offsets that satisfy
> the constraints (and no mixed models in the present case)
> 
> x <- runif(300,0,10)
> fitness <- rpois(300,10*dnorm(x,mean=5,sd=2))
> 
> objfn <- function(par) {
>    off <- (par[1]+par[2]*x)*x
>    fit <- glm(fitness ~ 1,family="poisson",offset=off)
>    - logLik(fit)
> }
> opt <- optim(c(1,-1),objfn,lower=c(0,-Inf),upper=c(Inf,0),method="L-BFGS-B")
> 
> mod <- glm(fitness ~
> 1,family="poisson",offset=with(opt,(par[1]+par[2]*x)*x))
> 
> plot(fitness ~ x)
> points(predict(mod,type="response") ~ x,col="red")
> 
> Best,
> 
> F.
> 
> > Is there a way to enforce such constraints in mixed models packages
> > available in R? I'm currently using lme4, but I'm happy to switch to any
> > other package.
> > 
> > Cheers,
> > Pierre.
> > 
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models