VGAM (version 1.1-9)

fill1: Creates a Matrix of Appropriate Dimension

Description

A support function for the argument xij, it generates a matrix of an appropriate dimension.

Usage

fill1(x, values = 0, ncolx = ncol(x))

Value

matrix(values, nrow=nrow(x), ncol=ncolx), i.e., a matrix consisting of values values, with the number of rows matching

x, and the default number of columns is the number of columns of x.

Arguments

x

A vector or matrix which is used to determine the dimension of the answer, in particular, the number of rows. After converting x to a matrix if necessary, the answer is a matrix of values values, of dimension nrow(x) by ncolx.

values

Numeric. The answer contains these values, which are recycled columnwise if necessary, i.e., as matrix(values, ..., byrow=TRUE).

ncolx

The number of columns of the returned matrix. The default is the number of columns of x.

Author

T. W. Yee

Details

The xij argument for vglm allows the user to input variables specific to each linear/additive predictor. For example, consider the bivariate logit model where the first/second linear/additive predictor is the logistic regression of the first/second binary response respectively. The third linear/additive predictor is log(OR) = eta3, where OR is the odds ratio. If one has ocular pressure as a covariate in this model then xij is required to handle the ocular pressure for each eye, since these will be different in general. [This contrasts with a variable such as age, the age of the person, which has a common value for both eyes.] In order to input these data into vglm one often finds that functions fill1, fill2, etc. are useful.

All terms in the xij and formula arguments in vglm must appear in the form2 argument too.

See Also

vglm.control, vglm, multinomial, Select.

Examples

fill1(runif(5))
fill1(runif(5), ncol = 3)
fill1(runif(5), val = 1, ncol = 3)

# Generate (independent) eyes data for the examples below; OR=1.
nn <- 1000  # Number of people
eyesdata <- data.frame(lop = round(runif(nn), 2),
                       rop = round(runif(nn), 2),
                       age = round(rnorm(nn, 40, 10)))
eyesdata <- transform(eyesdata,
  mop = (lop + rop) / 2,        # Mean ocular pressure
  op  = (lop + rop) / 2,        # Value unimportant unless plotting
# op  =  lop,                   # Choose this if plotting
  eta1 = 0 - 2*lop + 0.04*age,  # Linear predictor for left eye
  eta2 = 0 - 2*rop + 0.04*age)  # Linear predictor for right eye
eyesdata <- transform(eyesdata,
  leye = rbinom(nn, size=1, prob = logitlink(eta1, inverse = TRUE)),
  reye = rbinom(nn, size=1, prob = logitlink(eta2, inverse = TRUE)))

# Example 1. All effects are linear.
fit1 <- vglm(cbind(leye,reye) ~ op + age,
             family = binom2.or(exchangeable = TRUE, zero = 3),
             data = eyesdata, trace = TRUE,
             xij = list(op ~ lop + rop + fill1(lop)),
             form2 =  ~ op + lop + rop + fill1(lop) + age)
head(model.matrix(fit1, type = "lm"))   # LM model matrix
head(model.matrix(fit1, type = "vlm"))  # Big VLM model matrix
coef(fit1)
coef(fit1, matrix = TRUE)  # Unchanged with 'xij'
constraints(fit1)
max(abs(predict(fit1)-predict(fit1, new = eyesdata)))  # Okay
summary(fit1)
if (FALSE) {
plotvgam(fit1,
     se = TRUE)  # Wrong, e.g., coz it plots against op, not lop.
# So set op = lop in the above for a correct plot.
}

# Example 2. This uses regression splines on ocular pressure.
# It uses a trick to ensure common basis functions.
BS <- function(x, ...)
  sm.bs(c(x,...), df = 3)[1:length(x), , drop = FALSE]  # trick

fit2 <-
  vglm(cbind(leye,reye) ~ BS(lop,rop) + age,
       family = binom2.or(exchangeable = TRUE, zero = 3),
       data = eyesdata, trace = TRUE,
       xij = list(BS(lop,rop) ~ BS(lop,rop) +
                                BS(rop,lop) +
                                fill1(BS(lop,rop))),
       form2 = ~  BS(lop,rop) + BS(rop,lop) + fill1(BS(lop,rop)) +
                        lop + rop + age)
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)
summary(fit2)
[e-mail]
max(abs(predict(fit2) - predict(fit2, new = eyesdata)))  # Okay
predict(fit2, new = head(eyesdata))  # OR is 'scalar' as zero=3
max(abs(head(predict(fit2)) -
             predict(fit2, new = head(eyesdata))))  # Should be 0
if (FALSE) {
plotvgam(fit2, se = TRUE, xlab = "lop")  # Correct
}

# Example 3. Capture-recapture model with ephemeral and enduring
# memory effects. Similar to Yang and Chao (2005), Biometrics.
deermice <- transform(deermice, Lag1 = y1)
M.tbh.lag1 <-
  vglm(cbind(y1, y2, y3, y4, y5, y6) ~ sex + weight + Lag1,
       posbernoulli.tb(parallel.t = FALSE ~ 0,
                       parallel.b = FALSE ~ 0,
                       drop.b = FALSE ~ 1),
       xij = list(Lag1 ~ fill1(y1) + fill1(y2) + fill1(y3) +
                         fill1(y4) + fill1(y5) + fill1(y6) +
                         y1 + y2 + y3 + y4 + y5),
       form2 = ~ sex + weight + Lag1 +
                 fill1(y1) + fill1(y2) + fill1(y3) + fill1(y4) +
                 fill1(y5) + fill1(y6) +
                 y1 + y2 + y3 + y4 + y5 + y6,
       data = deermice, trace = TRUE)
coef(M.tbh.lag1)