mosaic.glm: Mosaic plots for fitted generalized linear and generalized...

In vcdExtra: 'vcd' Extensions and Additions

Mosaic plots for fitted generalized linear and generalized nonlinear models

Description

Produces mosaic plots (and other plots in the strucplot framework) for a log-linear model fitted with glm or for a generalized nonlinear model fitted with gnm.

These methods extend the range of strucplot visualizations well beyond the models that can be fit with loglm. They are intended for models for counts using the Poisson family (or quasi-poisson), but should be sensible as long as (a) the response variable is non-negative and (b) the predictors visualized in the strucplot are discrete factors.

Usage

## S3 method for class 'glm'
mosaic(x, formula = NULL, panel = mosaic, 
     type = c("observed", "expected"), 
     residuals = NULL, 
     residuals_type = c("pearson", "deviance", "rstandard"), 
     gp = shading_hcl, gp_args = list(), ...)
## S3 method for class 'glm'
sieve(x,  ...)
## S3 method for class 'glm'
assoc(x,  ...)

Arguments

x	A glm or gnm object. The response variable, typically a cell frequency, should be non-negative.
formula	A one-sided formula with the indexing factors of the plot separated by '+', determining the order in which the variables are used in the mosaic. A formula must be provided unless x$data inherits from class "table" – in which case the indexing factors of this table are used, or the factors in x$data (or model.frame(x) if x$data is an environment) exactly cross-classify the data – in which case this set of cross-classifying factors are used.
panel	Panel function used to draw the plot for visualizing the observed values, residuals and expected values. Currently, one of "mosaic", "assoc", or "sieve" in vcd.
type	A character string indicating whether the "observed" or the "expected" values of the table should be visualized by the area of the tiles or bars.
residuals	An optional array or vector of residuals corresponding to the cells in the data, for example, as calculated by residuals.glm(x), residuals.gnm(x).
residuals_type	If the residuals argument is NULL, residuals are calculated internally and used in the display. In this case, residual_type can be "pearson", "deviance" or "rstandard". Otherwise (when residuals is supplied), residuals_type is used as a label for the legend in the plot.
gp	Object of class "gpar", shading function or a corresponding generating function (see strucplot Details and shadings). Ignored if shade = FALSE.
gp_args	A list of arguments for the shading-generating function, if specified.
...	Other arguments passed to the panel function e.g., mosaic

Details

For both poisson family generalized linear models and loglinear models, standardized residuals provided by rstandard (sometimes called adjusted residuals) are often preferred because they have constant unit asymptotic variance.

The sieve and assoc methods are simple convenience interfaces to this plot method, setting the panel argument accordingly.

Value

The structable visualized by strucplot is returned invisibly.

Author(s)

Heather Turner, Michael Friendly, with help from Achim Zeileis

See Also

glm, gnm, plot.loglm, mosaic

Examples

GSStab <- xtabs(count ~ sex + party, data=GSS)
# using the data in table form
mod.glm1 <- glm(Freq ~ sex + party, family = poisson, data = GSStab)
res <- residuals(mod.glm1)    
std <- rstandard(mod.glm1)

# For mosaic.default(), need to re-shape residuals to conform to data
stdtab <- array(std, 
                dim=dim(GSStab), 
                dimnames=dimnames(GSStab))

mosaic(GSStab, 
       gp=shading_Friendly, 
       residuals=stdtab, 
       residuals_type="Std\nresiduals", 
       labeling = labeling_residuals)


# Using externally calculated residuals with the glm() object
mosaic.glm(mod.glm1, 
           residuals=std, 
           labeling = labeling_residuals, 
           shade=TRUE)

# Using residuals_type
mosaic.glm(mod.glm1, 
           residuals_type="rstandard", 
           labeling = labeling_residuals, shade=TRUE)

## Ordinal factors and structured associations
data(Mental)
xtabs(Freq ~ mental+ses, data=Mental)
long.labels <- list(set_varnames = c(mental="Mental Health Status", 
                                     ses="Parent SES"))

# fit independence model
# Residual deviance: 47.418 on 15 degrees of freedom
indep <- glm(Freq ~ mental+ses,
             family = poisson, data = Mental)

long.labels <- list(set_varnames = c(mental="Mental Health Status", 
                                     ses="Parent SES"))
mosaic(indep,
       residuals_type="rstandard", 
       labeling_args = long.labels, 
       labeling=labeling_residuals)

# or, show as a sieve diagram
mosaic(indep, 
       labeling_args = long.labels, 
       panel=sieve, 
       gp=shading_Friendly)

# fit linear x linear (uniform) association.  Use integer scores for rows/cols 
Cscore <- as.numeric(Mental$ses)
Rscore <- as.numeric(Mental$mental)

linlin <- glm(Freq ~ mental + ses + Rscore:Cscore,
                family = poisson, data = Mental)

mosaic(linlin,
       residuals_type="rstandard", 
       labeling_args = long.labels, 
       labeling=labeling_residuals, 
       suppress=1, 
       gp=shading_Friendly,
       main="Lin x Lin model")

##  Goodman Row-Column association model fits even better (deviance 3.57, df 8)
if (require(gnm)) {
Mental$mental <- C(Mental$mental, treatment)
Mental$ses <- C(Mental$ses, treatment)
RC1model <- gnm(Freq ~ ses + mental + Mult(ses, mental),
                family = poisson, data = Mental)

mosaic(RC1model,
       residuals_type="rstandard", 
       labeling_args = long.labels, 
       labeling=labeling_residuals, 
       suppress=1, 
       gp=shading_Friendly,
       main="RC1 model")
 }
 
 ############# UCB Admissions data, fit using glm()
 
structable(Dept ~ Admit+Gender,UCBAdmissions)
 
berkeley <- as.data.frame(UCBAdmissions)
berk.glm1 <- glm(Freq ~ Dept * (Gender+Admit), data=berkeley, family="poisson")
summary(berk.glm1)

mosaic(berk.glm1, 
       gp=shading_Friendly, 
       labeling=labeling_residuals, 
       formula=~Admit+Dept+Gender)

# the same, displaying studentized residuals; 
# note use of formula to reorder factors in the mosaic
mosaic(berk.glm1, 
       residuals_type="rstandard", 
       labeling=labeling_residuals, 
       shade=TRUE, 
	     formula=~Admit+Dept+Gender, 
	     main="Model: [DeptGender][DeptAdmit]")

## all two-way model
berk.glm2 <- glm(Freq ~ (Dept + Gender + Admit)^2, data=berkeley, family="poisson")
summary(berk.glm2)

mosaic.glm(berk.glm2, 
       residuals_type="rstandard", 
       labeling = labeling_residuals, 
       shade=TRUE,
	     formula=~Admit+Dept+Gender, 
	     main="Model: [DeptGender][DeptAdmit][AdmitGender]")

anova(berk.glm1, berk.glm2, test="Chisq")

# Add 1 df term for association of [GenderAdmit] only in Dept A
berkeley <- within(berkeley, 
                   dept1AG <- (Dept=='A')*(Gender=='Female')*(Admit=='Admitted'))
berkeley[1:6,]

berk.glm3 <- glm(Freq ~ Dept * (Gender+Admit) + dept1AG, data=berkeley, family="poisson")
summary(berk.glm3)

mosaic.glm(berk.glm3, 
           residuals_type="rstandard", 
           labeling = labeling_residuals, 
           shade=TRUE,
	         formula=~Admit+Dept+Gender, 
	         main="Model: [DeptGender][DeptAdmit] + DeptA*[GA]")

# compare models
anova(berk.glm1, berk.glm3, test="Chisq")
 

vcdExtra documentation built on Aug. 22, 2023, 9:11 a.m.