## housing.R Visualize models from example(housing, package="MASS")
# These examples fit a variety of models to the data(housing), giving a 4-way
# frequency table of 1681 individuals from the Copenhagen Housing Conditions
# Survey, classified by their Type of rental dwelling, perceived Influence on
# management of the property, and degree of Contact with other residents. The
# response variable here is Satisfaction of householders with their present
# housing circumstances.
library(vcdExtra)
library(MASS)
data(housing, package="MASS")
oldop <-options(contrasts = c("contr.treatment", "contr.poly"))
##########################
# Poisson models for Freq, equivalent to loglinear models
##########################
# Baseline model, with Satisfaction as a response
house.glm0 <- glm(Freq ~ Infl*Type*Cont + Sat, family = poisson,
data = housing)
modFit(house.glm0)
# labeling_args for mosaic()
largs <- list(set_varnames = c(
Infl="Influence on management",
Cont="Contact among residents",
Type="Type of dwelling",
Sat="Satisfaction"),
abbreviate=c(Type=3))
mosaic(house.glm0, labeling_args=largs, main='Baseline model: [ITC][Sat]')
# reorder variables in the mosaic, putting Sat last
mosaic(house.glm0, ~ Type+Infl+Cont+Sat, labeling_args=largs, main='Baseline model: [ITC][Sat]')
# what terms need to be added? Consider main effects, interactions of Sat with each other
MASS::addterm(house.glm0, ~. + Sat:(Infl+Type+Cont), test = "Chisq")
# add all two way terms with Satisfaction
house.glm1 <- update(house.glm0, . ~ . + Sat*(Infl+Type+Cont))
# did it get better?
anova(house.glm0, house.glm1, test="Chisq")
# plot it
mosaic(house.glm1, labeling_args=largs, main='Model [IS][TS][CS]', gp=shading_Friendly)
# Same model, fit by iterative proportional scaling
(house.loglm <- MASS::loglm(Freq ~ Infl*Type*Cont + Sat*(Infl+Type+Cont), data = housing))
# Can we drop any terms?
MASS::dropterm(house.glm1, test = "Chisq")
# Need to add any terms?
MASS::addterm(house.glm1, ~. + Sat:(Infl+Type+Cont)^2, test = "Chisq")
# add an interaction
house.glm2 <- update(house.glm1,
. ~ . + Sat:Infl:Type)
LRstats(house.glm0, house.glm1, house.glm2)
##########################
# Effect plots, for glm1 model
##########################
library(effects)
house.eff <-allEffects(house.glm1)
# show the interactions of Infl, Cont and Type with Sat
plot(house.eff, 'Infl:Sat', x.var='Sat', xlab="Satisfaction")
plot(house.eff, 'Infl:Sat', x.var='Infl', xlab="Influence")
# same plot in one panel, no std errors shown
plot(house.eff, 'Infl:Sat', x.var='Sat', xlab="Satisfaction", multiline=TRUE)
plot(house.eff, 'Cont:Sat', x.var='Sat', xlab="Satisfaction")
plot(house.eff, 'Type:Sat', x.var='Sat', xlab="Satisfaction")
##########################
# multinomial model
##########################
library(nnet)
# multinomial model, similar in spirit to house.glm1
(house.mult<- multinom(Sat ~ Infl + Type + Cont, weights = Freq,
data = housing))
# Do we need a more complex model?
house.mult2 <- multinom(Sat ~ Infl*Type*Cont, weights = Freq,
data = housing)
anova(house.mult, house.mult2)
# effect plots for multinomial model
house.effm <- allEffects(house.mult)
plot(house.effm, 'Infl', xlab='Influence on management', style="stacked",
main="Multinomial: Infl effect plot")
plot(house.effm, 'Cont', xlab='Contact among residents', style="stacked",
main="Multinomial: Cont effect plot")
plot(house.effm, 'Type', xlab='Type of dwelling', style="stacked",
main="Multinomial: Type effect plot")
##########################
# proportional odds model
##########################
(house.plr <- polr(Sat ~ Infl + Type + Cont,
data = housing, weights = Freq))
# Test proportional odds assumption by likelihood ratio test
# NB: multinom() objects do not have a df.residual component, so we have
# to use the difference in edf to get df for the test
pchisq(deviance(house.plr) - deviance(house.mult),
df = house.mult$edf -house.plr$edf, lower.tail = FALSE)
# try more complex models
house.plr2 <- stepAIC(house.plr, ~.^2)
house.plr2$anova
house.effp <- allEffects(house.plr)
plot(house.effp, 'Infl', xlab='Influence on management', style="stacked",
main="Proportional odds: Infl effect plot")
plot(house.effp, 'Cont', xlab='Contact among residents', style="stacked",
main="Proportional odds: Cont effect plot")
plot(house.effp, 'Type', xlab='Type of dwelling', style="stacked",
main="Proportional odds: Type effect plot")
options(oldop)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.