# Simulation Example
X1 <- rep(1,500)
X2 <- log(runif(500,0,30))
X3 <- log(runif(500,0,15))
X4 <- log(runif(500,10,20))
mui <- exp(-5 + 0.2*X2 - 0.03*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=gammahetero1(formula.mean,formula.shape)
a
Run the code above in your browser using DataLab