[R-sig-ME] GLMM with bounded parameters

[Ben Bolker]

  Late on this, but future reference: for glmer (not lmer), parameters
can be bounded by using the modular framework: following the example
from ?modular:

glmod <- glFormula(incidence/size ~ period + (1 | herd),
     data = cbpp, family = binomial, weight=size)
devfun <- do.call(mkGlmerDevfun, glmod)
opt <- optimizeGlmer(devfun)

At this point we can't quite use updateGlmerDevfun exactly as written

rho <- environment(devfun)
rho$nAGQ <- nAGQ  ## specify
## this is the default setting; you can substitute any other lower
##   bounds you like here ...
rho$lower <- c(rho$lower, rep.int(-Inf, length(rho$pp$beta0)))
rho$lp0 <- rho$pp$linPred(1)
rho$dpars <- seq_along(rho$pp$theta)
rho$baseOffset <- forceCopy(rho$resp$offset)
rho$GQmat <- GHrule(nAGQ)
rho$fac <- reTrms$flist[[1]]
devfun <- mkdevfun(rho, nAGQ)

  Then you're ready to finish:

opt <- optimizeGlmer(devfun, stage=2)
mkMerMod(environment(devfun), opt, glmod$reTrms, fr = glmod$fr)


On 16-12-19 10:44 AM, Pierre de Villemereuil wrote:
> 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
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>