
Last chance! 50% off unlimited learning
Sale ends in
grc(y, Rank = 1, Index.corner = 2:(1 + Rank),
str0 = 1, summary.arg = FALSE, h.step = 1e-04, ...)
rcim(y, family = poissonff, Rank = 0, Musual = NULL,
weights = NULL, which.linpred = 1,
Index.corner = ifelse(is.null(str0), 0, max(str0)) + 1:Rank,
rprefix = "Row.", cprefix = "Col.", iprefix = "X2.",
offset = 0, str0 = if (Rank) 1 else NULL,
summary.arg = FALSE, h.step = 0.0001,
rbaseline = 1, cbaseline = 1,
has.intercept = TRUE,
M = NULL,
rindex = 2:nrow(y),
cindex = 2:ncol(y),
iindex = 2:nrow(y),
...)
grc
a matrix of counts.
For rcim
a general matrix response depending on family
.
Output from table()
is acceptable; it is converted into a matrix.
Note that y
should be at least 3 bRank
interaction term.
Not all family functions are suitable or make sense.
All other linear/additmin(nrow(y), ncol(y))
}.
This is the dimension of the fit in terms of the interaction.
For grc()
this argument must be positive.
A value of 0 means no interactions (i.e., main effects only);
1:Musual
.Rank
integers.
These are used to store the Rank
by Rank
identity matrix in the
A
matrix; corner constraints are used.Rank = 0
, else
an integer from the set {1,...,min(nrow(y), ncol(y))
},
specifying the row that is used as the structural zero.
Passed into rrvglm.control
if TRUE
, a summary is returned.
If TRUE
, y
may be the output (fitted
object) of grc()
.summary.rrvglm()
. Only used when summary.arg = TRUE
.rrvglm.control()
.family
function
for an ordinary (univariate) response.
Then the number of linear predictors of the rcim()
fit is
usually the number of columns of y
multiplcindex
means column index, and these point to the columns
of y
which are part of the vector of linear/additive predictor
rindex
means row index, and these are similar to cindex
.
iindex
means interaction index, and these are similar to cindex
."grc"
, which currently is the same as
an "rrvglm"
object.
Currently,
a rank-0 rcim()
object is of class rcim0-class
,
else of class "rcim"
(this may change in the future).rcim()
is experimental at this stage and
may have bugs.
Quite a lot of expertise is needed when fitting and in its
interpretion thereof. For example, the constraint
matrices applies the reduced-rank regression to the first
(see which.linpred
)
linear predictor and the other linear predictors are intercept-only and
have a common value throughout the entire data set.
This means that, by default, family =
zipoissonff
is
appropriate but not
family =
zipoisson
.
Else set family =
zipoisson
and which.linpred = 2
.
To understand what is going on, do examine the constraint
matrices of the fitted object, and reconcile this with
Equations (4.3) to (4.5) of Yee and Hastie (2003).
The functions temporarily create a permanent data frame
called .grc.df
or .rcim.df
, which used
to be needed by summary.rrvglm()
. Then these
data frames are deleted before exiting the function.
If an error occurs, then the data frames may be present
in the workspace.
A %*% t(C)
,
the product of two `thin' matrices.
Indeed, A
and C
have Rank
columns.
By default, the first column and row of the interaction matrix
A %*% t(C)
is chosen
to be structural zeros, because str0 = 1
.
This means the first row of A
are all zeros.
This function uses options()$contrasts
to set up the row and
column indicator variables.
In particular, Equation (4.5) of Yee and Hastie (2003) is used.
These are called Row.
and Col.
(by default) followed
by the row or column number.
The function rcim()
is more general than grc()
.
Its default is a no-interaction model of grc()
, i.e.,
rank-0 and a Poisson distribution. This means that each
row and column has a dummy variable associated with it.
The first row and column is baseline.
The power of rcim()
is that many family
argument.
For example,
uninormal
fits something in between a 2-way
ANOVA with and without interactions,
alaplace2
with Rank = 0
is something like
medpolish
.
Others include
zipoissonff
,
negbinomial
.
Hopefully one day all family
argument, although the
result may not have meaning.
Yee, T. W. and Hadi, A. F. (2013) Row-column interaction models In preparation.
Goodman, L. A. (1981) Association models and canonical correlation in the analysis of cross-classifications having ordered categories. Journal of the American Statistical Association, 76, 320--334.
Documentation accompanying the
rrvglm
,
rrvglm.control
,
rrvglm-class
,
summary.grc
,
moffset
,
Rcim
,
Qvar
,
plotrcim0
,
multinomial
,
alcoff
,
crashi
,
auuc
,
olym08
,
olym12
,
poissonff
.grc1 <- grc(auuc) # Undergraduate enrolments at Auckland University in 1990
fitted(grc1)
summary(grc1)
grc2 <- grc(auuc, Rank = 2, Index.corner = c(2, 5))
fitted(grc2)
summary(grc2)
model3 <- rcim(auuc, Rank = 1, fam = multinomial,
M = ncol(auuc)-1, cindex = 2:(ncol(auuc)-1), trace = TRUE)
fitted(model3)
summary(model3)
# 2012 Summer Olympic Games in London
top10 <- head(olym12, n = 10)
grc.oly1 <- with(top10, grc(cbind(gold, silver, bronze)))
round(fitted(grc.oly1))
round(resid(grc.oly1, type = "response"), digits = 1) # Response residuals
summary(grc.oly1)
Coef(grc.oly1)
# Roughly median polish
rcim0 <- rcim(auuc, fam = alaplace2(tau = 0.5, intparloc = TRUE), trace = TRUE)
round(fitted(rcim0), digits = 0)
round(100 * (fitted(rcim0) - auuc) / auuc, digits = 0) # Discrepancy
depvar(rcim0)
round(coef(rcim0, matrix = TRUE), digits = 2)
Coef(rcim0, matrix = TRUE)
# constraints(rcim0)
names(constraints(rcim0))
# Compare with medpolish():
(med.a <- medpolish(auuc))
fv <- med.a$overall + outer(med.a$row, med.a$col, "+")
round(100 * (fitted(rcim0) - fv) / fv) # Hopefully should be all 0s
Run the code above in your browser using DataLab