# NOT RUN {
# 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)
# }
Run the code above in your browser using DataLab