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. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.