| fill1 {VGAM} | R Documentation |
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))
Arguments
x |
A vector or matrix which is used to determine
the dimension of the
answer, in particular, the number of rows.
After converting |
values |
Numeric.
The answer contains these values,
which are recycled columnwise if necessary, i.e.,
as |
ncolx |
The number of columns of the returned matrix.
The default is the number of columns of |
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.
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.
Note
The effect of the xij argument is after
other arguments such as
exchangeable and zero.
Hence xij does not affect constraint matrices.
Additionally, there are currently 3 other
identical fill1
functions, called fill2, fill3 and fill4;
if you need more then assign fill5 = fill6 = fill1 etc.
The reason for this is that if more than
one fill1 function is
needed then they must be unique.
For example, if M=4 then
xij = list(op ~ lop + rop + fill1(mop) + fill1(mop))
would reduce to
xij = list(op ~ lop + rop + fill1(mop)), whereas
xij = list(op ~ lop + rop + fill1(mop) + fill2(mop))
would retain
all M terms, which is needed.
In Examples 1 to 3 below, the xij argument
illustrates covariates
that are specific to a linear predictor.
Here, lop/rop are
the ocular pressures of the left/right eye
in an artificial dataset,
and mop is their mean.
Variables leye and reye
might be the presence/absence of a particular
disease on the LHS/RHS
eye respectively.
In Example 3, the xij argument illustrates fitting the
(exchangeable) model where there is a common smooth function
of the
ocular pressure. One should use regression splines since
s in vgam does not handle
the xij
argument. However, regression splines such as
bs and ns need
to have
the same basis functions here for both functions, and Example 3
illustrates a trick involving a function BS to obtain this,
e.g., same knots. Although regression splines create more than a
single column per term in the model matrix,
fill1(BS(lop,rop))
creates the required (same) number of columns.
Author(s)
T. W. Yee
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)
## Not run:
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.
## End(Not run)
# 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)
fit2@smart.prediction
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
## Not run:
plotvgam(fit2, se = TRUE, xlab = "lop") # Correct
## End(Not run)
# 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)