getContrasts | R Documentation |
Computes contrasts or scaled contrasts for a set of (non-eliminated)
parameters from a gnm
model, and computes standard errors
for the estimated contrasts. Where possible, quasi standard errors are
also computed.
getContrasts(model, set = NULL, ref = "first", scaleRef = "mean",
scaleWeights = NULL, dispersion = NULL, check = TRUE, ...)
model |
a model object of class |
set |
a vector of indices (numeric) or coefficient names
(character). If |
ref |
either a single numeric index, or a vector
of real numbers which sum to 1, or one of the character
strings |
scaleRef |
as for |
scaleWeights |
either |
dispersion |
either |
check |
|
... |
arguments to pass to other functions. |
The indices in set
must all be in 1:length(coef(object))
. If
set = NULL
, a dialog is presented for the selection
of indices (model coefficients).
For the set of coefficients selected, contrasts and their standard
errors are computed. A check is performed first on the estimability of
all such contrasts (if check = TRUE
) or on a specified subset
(if check
is a numeric index vector). The specific
contrasts to be computed are controlled by the choice of ref
:
this may be "first"
(the default), for contrasts with the first
of the selected coefficients, or "last"
for contrasts with the
last, or "mean"
for contrasts with the arithmetic mean of the
coefficients in the selected set; or it may be an arbitrary vector of
weights (summing to 1, not necessarily all non-negative) which specify
a weighted mean against which contrasts are taken; or it may be a
single index specifying one of the coefficients with which all
contrasts should be taken. Thus, for example, ref = 1
is
equivalent to ref = "first"
, and ref = c(1/3, 1/3, 1/3)
is equivalent to ref = "mean"
when there are three coefficients
in the selected set
.
The contrasts may be scaled by
\frac{1}{\sqrt{\sum_r v_r * d_r^2}}
where d_r
is a contrast of the r'th coefficient in set
with the reference level specified by scaleRef
and v
is a
vector of weights (of the same length as set
)
specified by scaleWeights
. If
scaleWeights
is NULL
(the default), scaleRef
is ignored and no scaling is performed. Other options for
scaleWeights
are "unit"
for weights equal to one and
"setLength"
for weights equal to the reciprocal of
length(set)
. If scaleRef
is the same as ref
, these options constrain the sum of squared
contrasts to 1 and length(set)
respectively.
Quasi-variances (and corresponding quasi standard errors) are
reported for unscaled contrasts where possible. These statistics are
invariant to the choice of ref
, see Firth (2003) or Firth and
Menezes (2004) for more details.
An object of class
qv
— see qvcalc
.
David Firth and Heather Turner
Firth, D (2003). Overcoming the reference category problem in the presentation of statistical models. Sociological Methodology 33, 1–18.
Firth, D and Menezes, R X de (2004). Quasi-variances. Biometrika 91, 65–80.
gnm
, se.gnm
,
checkEstimable
, qvcalc
,
ofInterest
### Unscaled contrasts ###
set.seed(1)
## Fit the "UNIDIFF" mobility model across education levels -- see ?yaish
unidiff <- gnm(Freq ~ educ*orig + educ*dest +
Mult(Exp(educ), orig:dest),
ofInterest = "[.]educ", family = poisson,
data = yaish, subset = (dest != 7))
## Examine the education multipliers (differences on the log scale):
unidiffContrasts <- getContrasts(unidiff, ofInterest(unidiff))
plot(unidiffContrasts,
main = "Unidiff multipliers (log scale): intervals based on
quasi standard errors",
xlab = "Education level", levelNames = 1:5)
### Scaled contrasts (elliptical contrasts) ###
set.seed(1)
## Goodman Row-Column association model fits well (deviance 3.57, df 8)
mentalHealth$MHS <- C(mentalHealth$MHS, treatment)
mentalHealth$SES <- C(mentalHealth$SES, treatment)
RC1model <- gnm(count ~ SES + MHS + Mult(SES, MHS),
family = poisson, data = mentalHealth)
## Row scores and column scores are both unnormalized in this
## parameterization of the model
## The scores can be normalized as in Agresti's eqn (9.15):
rowProbs <- with(mentalHealth, tapply(count, SES, sum) / sum(count))
colProbs <- with(mentalHealth, tapply(count, MHS, sum) / sum(count))
mu <- getContrasts(RC1model, pickCoef(RC1model, "[.]SES"),
ref = rowProbs, scaleRef = rowProbs,
scaleWeights = rowProbs)
nu <- getContrasts(RC1model, pickCoef(RC1model, "[.]MHS"),
ref = colProbs, scaleRef = colProbs,
scaleWeights = colProbs)
all.equal(sum(mu$qv[,1] * rowProbs), 0)
all.equal(sum(nu$qv[,1] * colProbs), 0)
all.equal(sum(mu$qv[,1]^2 * rowProbs), 1)
all.equal(sum(nu$qv[,1]^2 * colProbs), 1)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.