4

Does anyone know how to obtain the fitted phi parameter from the model output (mgcv GAM with a beta distribution)?

Referring to the betar gam example provided:

library(mgcv)
## Simulate some beta data...
set.seed(3);n<-400
dat <- gamSim(1,n=n)
mu <- binomial()$linkinv(dat$f/4-2)
phi <- .5
a <- mu*phi;b <- phi - a;
dat$y <- rbeta(n,a,b) 

bm <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=betar(link="logit"),data=dat)

The output shows a phi estimate under "Family:Beta regression"

> bm

Family: Beta regression(0.491) 
Link function: logit 

Formula:
y ~ s(x0) + s(x1) + s(x2) + s(x3)

Estimated degrees of freedom:
1.73 1.63 5.62 1.00  total = 10.98
StupidWolf
  • 45,075
  • 17
  • 40
  • 72
lmyt
  • 71
  • 2
  • Hi Luis, I can see this is your first time posting. Make sure to provide reproducible examples, instead of just referring to a webpage. I edited it for you this time as I am slightly familiar with this. Next time, please take note and write a more informed question :) – StupidWolf Nov 25 '19 at 23:54
  • Can you add some of your code of how you are trying to achieve this? – sarmahdi Nov 26 '19 at 01:19

1 Answers1

4

During the fitting, the gam function calls betar() which then estimates the phi parameter called theta. I quoted the description of the function below:

Description:

 Family for use with ‘gam’ or ‘bam’, implementing regression for
 beta distributed data on (0,1). A linear predictor controls the
 mean, mu of the beta distribution, while the variance is then
 mu(1-mu)/(1+phi), with parameter phi being estimated during
 fitting, alongside the smoothing parameters.

Usage:

 betar(theta = NULL, link = "logit",eps=.Machine$double.eps*100)
  Arguments:

 theta: the extra parameter (phi above).

Using the example you have, you just need to look under the gam object:

bm <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),family=betar(link="logit"),data=dat)
exp(bm$family$getTheta())
[1] 0.4913482
StupidWolf
  • 45,075
  • 17
  • 40
  • 72
  • I suppose this is true in general and that extra fitted parameters are listed in this 'Theta'. How does one know how the parameters are transformed (exp in the case above)? I could not find this in the mgcv documentation. – lmyt Nov 26 '19 at 15:50
  • Hey @LuisMier-y-Teran, I looked under the source code for the family, betar, https://github.com/cran/mgcv/blob/master/R/efam.r, and you can see that most of the time, before they feed it into the beta-binomial distribution, the exponential is taken. – StupidWolf Nov 26 '19 at 18:18
  • 3
    +1 if you look at the `getTheta` function, you can do `getTheta(trans = TRUE)` to achieve the same thing as explicitly `exp()`-ing the result of the call to `getTheta()`. – Gavin Simpson Nov 26 '19 at 20:41
  • 1
    Yes @ReinstateMonica-G.Simpson, you are right. Thanks for pointing it out. – StupidWolf Nov 26 '19 at 21:07