### Modelling the mean in terms of x2, two responses.
set.seed(2017022101)
nn <- 80
x2 <- runif(nn)
mu <- exp(2 + 0.5 * x2)
# Shape and rate parameters in terms of 'mu'
shape <- rep(exp(1), nn)
rate <- gammaRMlink(theta = log(mu), shape = shape,
inverse = TRUE, deriv = 0)
# Generating some random data
y1 <- rgamma(n = nn, shape = shape, rate = rate)
gdata <- data.frame(x2 = x2, y1 = y1)
rm(y1)
# lmu = "gammaRMlink" replaces lshape, whilst lrate = "loglink"
fit1 <- vglm(cbind(y1, y1) ~ x2,
gammaRff(lmu = "gammaRMlink", lss = TRUE, zero = "shape"),
data = gdata, trace = TRUE, crit = "log")
coef(fit1, matrix = TRUE)
summary(fit1)
# Comparing fitted values with true values.
compare1 <- cbind(fitted.values(fit1)[, 1, drop = FALSE], mu)
colnames(compare1) <- c("Fitted.vM1", "mu")
head(compare1)
### Mimicking gammaR. Note that lmu = NULL.
fit2 <- vglm(y1 ~ x2, gammaRff(lmu = NULL, lrate = "loglink",
lshape = "loglink", lss = FALSE, zero = "shape"),
data = gdata, trace = TRUE, crit = "log")
# Compare fitted values with true values.
compare2 <- with(gdata, cbind(fitted.values(fit2), y1, mu))
colnames(compare2) <- c("Fitted.vM2", "y", "mu")
head(compare2)
### Fitted values -- Model1 vs Fitted values -- Model2
fit1vsfit2 <- cbind(fitted.values(fit1)[, 1, drop = FALSE],
fitted.values(fit2))
colnames(fit1vsfit2) <- c("Fitted.vM1", "Fitted.vM2")
head(fit1vsfit2)
### Use gamma2()
fit3 <- vglm(y1 ~ x2, gamma2,
data = gdata, trace = TRUE, crit = "log")
fit1.fit3 <- cbind(fitted.values(fit1)[, 1, drop = FALSE],
fitted.values(fit2), fitted.values(fit3))
colnames(fit1.fit3) <- c("Fitted.vM1", "Fitted.vM2", "Fitted.vM3")
head(fit1.fit3)
Run the code above in your browser using DataLab