glm() scale parameters and predicted Values

Peter Maclean <pmaclean2011 <at> yahoo.com> writes:

In glm() you can use the summary() function to recover
the shape parameter (the reciprocal of the
dispersion parameter). How do you recover the scale parameter? 
Also, in the given example, how I estimate
and save the geometric mean of the predicted values? 
For a simple model you can use fitted() or predicted()
functions. I will appreciate any help. 
?
?
?
#Call required R packages
require(plyr)? 
require(stats) 
require(fitdistrplus)
require(MASS)
#Grouped vector
n <- c(1:10)
yr <-c(1:10)
ny <- list(yr=yr,n=n)
require(utils)
ny <- expand.grid(ny) 
y = rgamma(100, shape=1.5, rate = 1, scale = 2)

It's a bit odd to specify both the rate and scale parameters,
which are redundant.  I would have guessed that the
rate parameter would dominate, but it looks like the
scale parameter dominates:

set.seed(1001); rgamma(1,shape=1,rate=2)

[1] 1.622577

set.seed(1001); rgamma(1,shape=1,scale=2)

[1] 6.490306

set.seed(1001); rgamma(1,shape=1,rate=2,scale=2)

[1] 6.490306

set.seed(1001); rgamma(1,shape=1,scale=2,rate=2)

[1] 6.490306

  (I know, I could go look at the source, but it's
a .Internal() function and I'm feeling lazy ...)

Gdata <- cbind(ny,y)
Gdata2<- Gdata
Gdata$x1 <- cos((3.14*yr)/365.25) 
Gdata$x2 <- sin((3.14*yr)/365.25) 
#Fitting Generalized Linear Models 
Gdata <- split(Gdata,Gdata$n)
FGLM <- lapply(Gdata, function(x){
????????????? m <- as.numeric(x$y)
????????????? x1 <- m <- as.numeric(x$x1)
????????????? x2 <- m <- as.numeric(x$x2)
????????????? summary(glm(m~1+x1+x2, family=Gamma),dispersion=NULL) 
?????????????? })

Note that you have simulated a null model (the data don't
depend on the covariates).

#Save the results of the estimated parameters
str(FGLM,no.list = TRUE)
SFGLMC<- ldply(FGLM, function(x) x$coefficients)
SFGLMD<- ldply(FGLM, function(x) x$dispersion)
GLM <- cbind(SFGLMC,SFGLMD)

Which scale parameter do you mean?  In a real
model it will vary with x1 and x2.  Let's try an
experiment first:

set.seed(1001)
z <- rgamma(10000,shape=2,scale=3)
g1 <- glm(z~1,family=Gamma)
summary(g1)

Call:
glm(formula = z ~ 1, family = Gamma)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8434  -0.6579  -0.1702   0.3081   2.6220  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.167013   0.001189   140.5   <2e-16 ***
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 

(Dispersion parameter for Gamma family taken to be 0.5066566)

    Null deviance: 5419.8  on 9999  degrees of freedom
Residual deviance: 5419.8  on 9999  degrees of freedom
AIC: 53526

  Here the intercept estimate is 0.167 , which is very
nearly 1/6 = 1/(shape*scale) -- i.e. the Gamma GLM is
parameterized in terms of the mean (on the inverse
scale).  If you want to recover
the scale parameter for the intercept case, then

summary(g1)$dispersion/coef(g1)[1]

should be pretty good.

As for the geometric means -- that's pretty tricky.
*If* you used a log link, then the predicted values on
the link scale (i.e. predict(g1,type="link")) would indeed
give you the geometric means.  On the inverse scale, though,
you would have to do a bit of finagling.