# Example 1.
pneumo <- transform(pneumo, let = log(exposure.time))
vglm(cbind(normal, mild, severe) ~ let, multinomial, pneumo,
crit = "coef", step = 0.5, trace = TRUE, epsil = 1e-8,
maxit = 40)
# Example 2. The use of the xij argument (simple case).
ymat <- rdiric(n <- 1000, shape = rep(exp(2), len = 4))
mydat <- data.frame(x1 = runif(n), x2 = runif(n), x3 = runif(n),
x4 = runif(n),
z1 = runif(n), z2 = runif(n), z3 = runif(n),
z4 = runif(n))
mydat <- transform(mydat, X = x1, Z = z1)
mydat <- round(mydat, digits = 2)
fit2 <- vglm(ymat ~ X + Z,
dirichlet(parallel = TRUE), mydat, trace = TRUE,
xij = list(Z ~ z1 + z2 + z3 + z4,
X ~ x1 + x2 + x3 + x4),
form2 = ~ Z + z1 + z2 + z3 + z4 +
X + x1 + x2 + x3 + x4)
head(model.matrix(fit2, type = "lm")) # LM model matrix
head(model.matrix(fit2, type = "vlm")) # Big VLM model matrix
coef(fit2)
coef(fit2, matrix = TRUE)
max(abs(predict(fit2)-predict(fit2, new = mydat))) # Predicts ok
summary(fit2)
if (FALSE) {
# plotvgam(fit2, se = TRUE, xlab = "x1", which.term = 1) # Bug!
# plotvgam(fit2, se = TRUE, xlab = "z1", which.term = 2) # Bug!
plotvgam(fit2, xlab = "x1") # Correct
plotvgam(fit2, xlab = "z1") # Correct
}
# Example 3. The use of the xij argument (complex case).
set.seed(123)
coalminers <-
transform(coalminers,
Age = (age - 42) / 5,
dum1 = round(runif(nrow(coalminers)), digits = 2),
dum2 = round(runif(nrow(coalminers)), digits = 2),
dum3 = round(runif(nrow(coalminers)), digits = 2),
dumm = round(runif(nrow(coalminers)), digits = 2))
BS <- function(x, ..., df = 3)
sm.bs(c(x,...), df = df)[1:length(x),,drop = FALSE]
NS <- function(x, ..., df = 3)
sm.ns(c(x,...), df = df)[1:length(x),,drop = FALSE]
# Equivalently...
BS <- function(x, ..., df = 3)
head(sm.bs(c(x,...), df = df), length(x), drop = FALSE)
NS <- function(x, ..., df = 3)
head(sm.ns(c(x,...), df = df), length(x), drop = FALSE)
fit3 <- vglm(cbind(nBnW,nBW,BnW,BW) ~ Age + NS(dum1, dum2),
fam = binom2.or(exchangeable = TRUE, zero = 3),
xij = list(NS(dum1, dum2) ~ NS(dum1, dum2) +
NS(dum2, dum1) +
fill1(NS( dum1))),
form2 = ~ NS(dum1, dum2) + NS(dum2, dum1) +
fill1(NS(dum1)) +
dum1 + dum2 + dum3 + Age + age + dumm,
data = coalminers, trace = TRUE)
head(model.matrix(fit3, type = "lm")) # LM model matrix
head(model.matrix(fit3, type = "vlm")) # Big VLM model matrix
coef(fit3)
coef(fit3, matrix = TRUE)
if (FALSE) {
plotvgam(fit3, se = TRUE, lcol = 2, scol = 4, xlab = "dum1")
}
Run the code above in your browser using DataLab