# Find more examples and instructions in the package vignette:
# vignette('hglm')
require(hglm)
# --------------------- #
# semiconductor example #
# --------------------- #
data(semiconductor)
m11 <- hglm(fixed = y ~ x1 + x3 + x5 + x6,
random = ~ 1|Device,
family = Gamma(link = log),
disp = ~ x2 + x3, data = semiconductor)
summary(m11)
plot(m11, cex = .6, pch = 1,
cex.axis = 1/.6, cex.lab = 1/.6,
cex.main = 1/.6, mar = c(3, 4.5, 0, 1.5))
# ------------------- #
# redo it using hglm2 #
# ------------------- #
m12 <- hglm2(y ~ x1 + x3 + x5 + x6 + (1|Device),
family = Gamma(link = log),
disp = ~ x2 + x3, data = semiconductor)
summary(m12)
# -------------------------- #
# redo it using matrix input #
# -------------------------- #
attach(semiconductor)
m13 <- hglm(y = y, X = model.matrix(~ x1 + x3 + x5 + x6),
Z = kronecker(diag(16), rep(1, 4)),
X.disp = model.matrix(~ x2 + x3),
family = Gamma(link = log))
summary(m13)
# --------------------- #
# verbose & likelihoods #
# --------------------- #
m14 <- hglm(fixed = y ~ x1 + x3 + x5 + x6,
random = ~ 1|Device,
family = Gamma(link = log),
disp = ~ x2 + x3, data = semiconductor,
verbose = TRUE, calc.like = TRUE)
summary(m14)
# --------------------------------------------- #
# simulated example with 2 random effects terms #
# --------------------------------------------- #
## Not run:
# set.seed(911)
# x1 <- rnorm(100)
# x2 <- rnorm(100)
# x3 <- rnorm(100)
# z1 <- factor(rep(LETTERS[1:10], rep(10, 10)))
# z2 <- factor(rep(letters[1:5], rep(20, 5)))
# Z1 <- model.matrix(~ 0 + z1)
# Z2 <- model.matrix(~ 0 + z2)
# u1 <- rnorm(10, 0, sqrt(2))
# u2 <- rnorm(5, 0, sqrt(3))
# y <- 1 + 2*x1 + 3*x2 + Z1%*%u1 + Z2%*%u2 + rnorm(100, 0, sqrt(exp(x3)))
# dd <- data.frame(x1 = x1, x2 = x2, x3 = x3, z1 = z1, z2 = z2, y = y)
#
# (m21 <- hglm(X = cbind(rep(1, 100), x1, x2), y = y, Z = cbind(Z1, Z2),
# RandC = c(10, 5)))
# summary(m21)
# plot(m21)
#
# # m21 is the same as:
# (m21b <- hglm(X = cbind(rep(1, 100), x1, x2), y = y, Z = cbind(Z1, Z2),
# rand.family = list(gaussian(), gaussian()), RandC = c(10, 5))
#
# (m22 <- hglm2(y ~ x1 + x2 + (1|z1) + (1|z2), data = dd, vcovmat = TRUE))
# image(m22$vcov, main = 'Variance-covariance Matrix')
# summary(m22)
# plot(m22)
#
# m31 <- hglm2(y ~ x1 + x2 + (1|z1) + (1|z2), disp = ~ x3, data = dd)
# print (m31)
# summary(m31)
# plot(m31)
#
# # ------------------------------- #
# # Markov random field (MRF) model #
# # ------------------------------- #
# data(cancer)
# logE <- log(E)
# X11 <- model.matrix(~Paff)
# m41 <- hglm(X = X11, y = O, Z = diag(length(O)),
# family = poisson(), rand.family = CAR(D = nbr),
# offset = logE, conv = 1e-9, maxit = 200, fix.disp = 1)
# summary(m41)
#
# data(ohio)
# m42 <- hglm(fixed = MedianScore ~ 1,
# random = ~ 1 | district,
# rand.family = CAR(D = ohioDistrictDistMat),
# data = ohioMedian)
# summary(m42)
# require(sp)
# districtShape <- as.numeric(substr(as.character(ohioShape@data$UNSDIDFP), 3, 7))
# CARfit <- matrix(m42$ranef + m42$fixef, dimnames = list(rownames(ohioDistrictDistMat), NULL))
# ohioShape@data$CAR <- CARfit[as.character(districtShape),]
# ohioShape@data$CAR[353] <- NA # remove estimate of Lake Erie
# spplot(ohioShape, zcol = "CAR", main = "Fitted values from CAR",
# col.regions = heat.colors(1000)[1000:1], cuts = 1000)
# ## End(Not run)
Run the code above in your browser using DataLab