Last data update: 2014.03.03
R: Gammareg
Gammareg
Description
Function to do Classic Gamma Regression: joint mean and shape modeling
Usage
Gammareg(formula1,formula2,meanlink)
Arguments
formula1
object of class matrix, with the dependent variable.
formula2
object of class matrix, with the variables for modelling the mean.
meanlink
links for the mean. The default links is the link log. The link identity is also allowed as admisible value.
Details
The classic gamma regression allow the joint modelling of mean and shape parameters of a gamma distributed variable,
as is proposed in Cepeda (2001), using the Fisher Socring algorithm, with two differentes link for the mean: log and identity, and log link for the shape.
Value
object of class bayesbetareg with:
coefficients
object of class matrix with the estimated coefficients of beta and gamma.
desvB
object of class matrix with the estimated covariances of beta.
desvG
object of class matrix with the estimated covariances of gamma.
interv
object of class matrix with the estimated confidence intervals of beta and gamma.
AIC
the AIC criteria.
iteration
numbers of iterations to convergence.
convergence
value of convergence obtained.
call
Call.
Author(s)
Martha Corrales martha.corrales@usa.edu.co
Edilberto Cepeda-Cuervo ecepedac@unal.edu.co
References
1. Cepeda-Cuervo, E. (2001). Modelagem da variabilidade em modelos lineares generalizados. Unpublished Ph.D. tesis. Instituto de Matem<c3><a1>ticas.
Universidade Federal do R<c3><ad>o do Janeiro.
//http://www.docentes.unal.edu.co/ecepedac/docs/MODELAGEM20DA20VARIABILIDADE.pdf.
http://www.bdigital.unal.edu.co/9394/.
2. McCullagh, P. and Nelder, N.A. (1989). Generalized Linear Models. Second Edition. Chapman and Hall.
Examples
#
num.killed <- c(7,59,115,149,178,229,5,43,76,4,57,83,6,57,84)
size.sam <- c(1,2,3,3,3,3,rep(1,9))*100
insecticide <- c(4,5,8,10,15,20,2,5,10,2,5,10,2,5,10)
insecticide.2 <- insecticide^2
synergist <- c(rep(0,6),rep(3.9,3),rep(19.5,3),rep(39,3))
par(mfrow=c(2,2))
plot(density(num.killed/size.sam),main="")
boxplot(num.killed/size.sam)
plot(insecticide,num.killed/size.sam)
plot(synergist,num.killed/size.sam)
mean.for <- (num.killed/size.sam) ~ insecticide + insecticide.2
dis.for <- ~ synergist + insecticide
res=Gammareg(mean.for,dis.for,meanlink="ide")
summary(glm((num.killed/size.sam) ~ insecticide + insecticide.2,family=Gamma("log")))
summary(res)
# Simulation Example
X1 <- rep(1,500)
X2 <- runif(500,0,30)
X3 <- runif(500,0,15)
X4 <- runif(500,10,20)
mui <- 15 + 2*X2 + 3*X3
alphai <- exp(0.2 + 0.1*X2 + 0.3*X4)
Y <- rgamma(500,shape=alphai,scale=mui/alphai)
X <- cbind(X1,X2,X3)
Z <- cbind(X1,X2,X4)
formula.mean= Y~X2+X3
formula.shape= ~X2+X4
a=Gammareg(formula.mean,formula.shape,meanlink="ide")
summary(a)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(Gammareg)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/Gammareg/Gammareg.Rd_%03d_medium.png", width=480, height=480)
> ### Name: Gammareg
> ### Title: Gammareg
> ### Aliases: Gammareg
> ### Keywords: Classic estimation Fisher Scoring joint modeling Gamma
> ### regression Mean and shape parameters
>
> ### ** Examples
>
>
> #
>
> num.killed <- c(7,59,115,149,178,229,5,43,76,4,57,83,6,57,84)
> size.sam <- c(1,2,3,3,3,3,rep(1,9))*100
> insecticide <- c(4,5,8,10,15,20,2,5,10,2,5,10,2,5,10)
> insecticide.2 <- insecticide^2
> synergist <- c(rep(0,6),rep(3.9,3),rep(19.5,3),rep(39,3))
>
> par(mfrow=c(2,2))
> plot(density(num.killed/size.sam),main="")
> boxplot(num.killed/size.sam)
> plot(insecticide,num.killed/size.sam)
> plot(synergist,num.killed/size.sam)
>
>
> mean.for <- (num.killed/size.sam) ~ insecticide + insecticide.2
> dis.for <- ~ synergist + insecticide
>
> res=Gammareg(mean.for,dis.for,meanlink="ide")
>
> summary(glm((num.killed/size.sam) ~ insecticide + insecticide.2,family=Gamma("log")))
Call:
glm(formula = (num.killed/size.sam) ~ insecticide + insecticide.2,
family = Gamma("log"))
Deviance Residuals:
Min 1Q Median 3Q Max
-0.8439 -0.4926 -0.1526 0.2123 0.8703
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.48260 0.47761 -7.292 9.58e-06 ***
insecticide 0.52757 0.11376 4.638 0.000573 ***
insecticide.2 -0.01911 0.00546 -3.500 0.004381 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 0.3815866)
Null deviance: 11.9764 on 14 degrees of freedom
Residual deviance: 4.1655 on 12 degrees of freedom
AIC: -3.8819
Number of Fisher Scoring iterations: 8
> summary(res)
################################################################
### Classic Gamma Regression ###
################################################################
Call:
Gammareg(formula1 = mean.for, formula2 = dis.for, meanlink = "ide")
Estimate L.Intv U.Intv
beta.(Intercept) -0.129697 -0.261399 0.002
beta.insecticide 0.119233 0.068680 0.170
beta.insecticide.2 -0.004089 -0.006540 -0.002
gamma.(Intercept) 0.439037 0.182929 0.695
gamma.synergist 0.015735 0.010504 0.021
gamma.insecticide 0.162918 0.149205 0.177
Covariance Matrix for Beta:
(Intercept) insecticide insecticide.2
(Intercept) 3.580541e-03 -1.213089e-03 5.459089e-05
insecticide -1.213089e-03 5.275563e-04 -2.494403e-05
insecticide.2 5.459089e-05 -2.494403e-05 1.239960e-06
Covariance Matrix for Gamma:
(Intercept) synergist insecticide
(Intercept) 0.0138168350 -2.122943e-04 -7.266857e-04
synergist -0.0002122943 5.764239e-06 1.080063e-05
insecticide -0.0007266857 1.080063e-05 3.961207e-05
AIC:
[1] -577.8824
Iteration:
[1] 59
Convergence:
[1] 5.61794e-06
>
> # Simulation Example
>
> X1 <- rep(1,500)
> X2 <- runif(500,0,30)
> X3 <- runif(500,0,15)
> X4 <- runif(500,10,20)
> mui <- 15 + 2*X2 + 3*X3
> alphai <- exp(0.2 + 0.1*X2 + 0.3*X4)
> Y <- rgamma(500,shape=alphai,scale=mui/alphai)
> X <- cbind(X1,X2,X3)
> Z <- cbind(X1,X2,X4)
> formula.mean= Y~X2+X3
> formula.shape= ~X2+X4
> a=Gammareg(formula.mean,formula.shape,meanlink="ide")
> summary(a)
################################################################
### Classic Gamma Regression ###
################################################################
Call:
Gammareg(formula1 = formula.mean, formula2 = formula.shape, meanlink = "ide")
Estimate L.Intv U.Intv
beta.(Intercept) 15.3913 14.9044 15.878
beta.X2 2.0001 1.9763 2.024
beta.X3 2.9851 2.9399 3.030
gamma.(Intercept) 0.3195 0.3188 0.320
gamma.X2 0.1048 0.1048 0.105
gamma.X4 0.2866 0.2866 0.287
Covariance Matrix for Beta:
(Intercept) X2 X3
(Intercept) 0.061409236 -2.264497e-03 -2.365192e-03
X2 -0.002264497 1.460311e-04 -3.309234e-05
X3 -0.002365192 -3.309234e-05 5.280078e-04
Covariance Matrix for Gamma:
(Intercept) X2 X4
(Intercept) 1.311470e-07 -1.008708e-09 -5.763330e-09
X2 -1.008708e-09 3.772879e-11 4.617491e-12
X4 -5.763330e-09 4.617491e-12 3.084264e-10
AIC:
[1] 7113416
Iteration:
[1] 11
Convergence:
[1] 3.365878e-10
>
>
>
>
>
>
> dev.off()
null device
1
>