# [R] Bayesian functions for mle2 object

Ben Bolker bbolker at gmail.com
Thu Sep 1 03:02:31 CEST 2011

```On 11-09-05 02:13 PM, Billy wrote:
> Dear Dr. Ben Bolker and 'JLucke',
>
> For example, using the gls() function, I receive an error message and
> I'm not sure I'm writing the arguments in the right way.

Well, (if you still need help with this -- your comments below
seem like you're making progress with a different approach)
what did you try, and what is the error message?  We can't help
you without details ...

> Instead, I thought about my original models and realized that I was
> modelling variance as a linear function of the predictor variable, which
> could drop off to zero values. Changing
>
> sd = c0 + c1*x
>
> to
>
> sd = c0*x^c1
>
> avoids the zero values and all problematic models have converged. Paying
> attention on the estimates, they also make sense.

Good.

>
> The new problem now (and probably due to my weak Mathematical skills) is
> that one set o candidate models includes models that consider the effect
> of not only one, but two predictor variables on the response (y).
>
> How could be the right way to model that?
>
> sd = c0 * (x ^ c1) * (w ^ c2)
>
> or
>
> sd = (c0 * x ^ c1) + (c0 * w ^ c2) ?
>
> In which c0, c1, and c2 are constant parameters, and x and w are
> different predictor variables.
>
> Thanks again
>
> Billy

It's not clear that either of them is necessarily more "right" than
the others, but you could (much) more easily implement the first in
gls(); take a look at ?varComb ... you would use something like

weights=varComb(varPower(form=~x),varPower(form=~w))

>
> On Mon, Aug 29, 2011 at 3:50 PM, Ben Bolker <bbolker at gmail.com
> <mailto:bbolker at gmail.com>> wrote:
>
>     Billy.Requena <billy.requena <at> gmail.com <http://gmail.com>> writes:
>
>     >
>     > Hi everybody,
>     >
>     > I'm interested in evaluating the effect of a continuous variable
>     on the mean
>     > and/or the variance of my response variable. I have built functions
>     > expliciting these and used the 'mle2' function to estimate the
>     coefficients,
>     > as follows:
>     >
>     > func.1 <- function(m=62.9, c0=8.84, c1=-1.6)
>     >       {
>     >       s <- c0+c1*(x)
>     >       -sum(dnorm(y, mean=m, sd=s,log=T))
>     >       }
>     >
>     > m1 <- mle2(func.1, method="SANN")
>     >
>     > However, the estimation of the effect of x on the variance of y
>     usually has
>     > dealt some troubles, resulting in no convergencies or sd of estimates
>     > extremely huge. I tried using different optimizers, but I still
>     faced the
>     > some problems.
>     >
>     > When I had similar troubles in 'GLMM' statistical universe, I used
>     bayesian
>     > functions to solve this problem, enjoyning the flexibility of
>     different
>     > start points to reach the maximum likelihood estimates. However, I
>     have no
>     > idea which package or which function to use to solve the specific
>     problem
>     > I'm facing now.
>     > Does anyone have a clue?
>
>      Unless I'm missing something, you can fit this model
>     (more easily) in gls() from the nlme package, which allows models
>     for heteroscedasticity.  See ?nlme::varConstPower
>
>      gls(y~1,weights=varPower(power=1,form=~x),data)
>
>      This gives you a standard deviation proportional to (t1+|v|);
>     that is, if the baseline residual standard deviation is S, then
>     the standard deviation is S*(t1+|v|), so S would correspond to
>
>       Ben Bolker
>
>     ______________________________________________
>     R-help at r-project.org <mailto:R-help at r-project.org> mailing list
>     https://stat.ethz.ch/mailman/listinfo/r-help
>     http://www.R-project.org/posting-guide.html
>     and provide commented, minimal, self-contained, reproducible code.
>
>
>
>
> --
> Gustavo Requena
> PhD student - Laboratory of Arthropod Behavior and Evolution