> library(MASS) > data(housing) # Loads data set > head(housing) Sat Infl Type Cont Freq 1 Low Low Tower Low 21 2 Medium Low Tower Low 21 3 High Low Tower Low 28 4 Low Medium Tower Low 34 5 Medium Medium Tower Low 22 6 High Medium Tower Low 36 > ?housing starting httpd help server ... done > class(housing$Sat) # Treating as ordered factor [1] "ordered" "factor" > head(housing$Sat) [1] Low Medium High Low Medium High Levels: Low < Medium < High > housing$Sat <- factor(housing$Sat,ordered=FALSE) > class(housing$Infl) [1] "factor" > attach(housing) > sum(Freq) # Sample size [1] 1681 > House <- xtabs(Freq ~ Type+Infl+Cont+Sat, housing) Warning messages: 1: In Ops.factor(Type, Infl) : ‘+’ not meaningful for factors 2: In Ops.factor(Type + Infl, Cont) : ‘+’ not meaningful for factors 3: In Ops.factor(Type + Infl + Cont, Sat) : ‘+’ not meaningful for factors > House , , Cont = Low, Sat = Low Infl Type Low Medium High Tower 21 34 10 Apartment 61 43 26 Atrium 13 8 6 Terrace 18 15 7 , , Cont = High, Sat = Low Infl Type Low Medium High Tower 14 17 3 Apartment 78 48 15 Atrium 20 10 7 Terrace 57 31 5 , , Cont = Low, Sat = Medium Infl Type Low Medium High Tower 21 22 11 Apartment 23 35 18 Atrium 9 8 7 Terrace 6 13 5 , , Cont = High, Sat = Medium Infl Type Low Medium High Tower 19 23 5 Apartment 46 45 25 Atrium 23 22 10 Terrace 23 21 6 , , Cont = Low, Sat = High Infl Type Low Medium High Tower 28 36 36 Apartment 17 40 54 Atrium 10 12 9 Terrace 7 13 11 , , Cont = High, Sat = High Infl Type Low Medium High Tower 37 40 23 Apartment 43 86 62 Atrium 20 24 21 Terrace 13 13 13 > mosaicplot(House, color=TRUE) > House <- xtabs(Freq ~ Sat+Type+Infl+Cont, housing) > mosaicplot(House, type="deviance", shade=T) > tab1 <- tapply(Freq,list(Cont,Sat),sum) Warning message: In Ops.factor(Type, Infl) : ‘+’ not meaningful for factors > tab1 <- tapply(Freq,list(Cont,Sat),sum) > tab1 Low Medium High Low 262 178 273 High 305 268 395 > mod1 <- glm(Freq ~ Cont+Sat, family=poisson) > mod1 Call: glm(formula = Freq ~ Cont + Sat, family = poisson) Coefficients: (Intercept) ContHigh SatMedium SatHigh 2.9978 0.3058 -0.2400 0.1639 Degrees of Freedom: 71 Total (i.e. Null); 68 Residual Null Deviance: 833.7 Residual Deviance: 750.2 AIC: 1099 > mod1$fitted 1 2 3 4 5 6 20.04120 15.76433 23.61114 20.04120 15.76433 23.61114 7 8 9 10 11 12 20.04120 15.76433 23.61114 20.04120 15.76433 23.61114 13 14 15 16 17 18 20.04120 15.76433 23.61114 20.04120 15.76433 23.61114 19 20 21 22 23 24 20.04120 15.76433 23.61114 20.04120 15.76433 23.61114 25 26 27 28 29 30 20.04120 15.76433 23.61114 20.04120 15.76433 23.61114 31 32 33 34 35 36 20.04120 15.76433 23.61114 20.04120 15.76433 23.61114 37 38 39 40 41 42 27.20880 21.40234 32.05552 27.20880 21.40234 32.05552 43 44 45 46 47 48 27.20880 21.40234 32.05552 27.20880 21.40234 32.05552 49 50 51 52 53 54 27.20880 21.40234 32.05552 27.20880 21.40234 32.05552 55 56 57 58 59 60 27.20880 21.40234 32.05552 27.20880 21.40234 32.05552 61 62 63 64 65 66 27.20880 21.40234 32.05552 27.20880 21.40234 32.05552 67 68 69 70 71 72 27.20880 21.40234 32.05552 27.20880 21.40234 32.05552 > mod1.fits <- tapply(mod1$fitted, list(Cont,Sat), sum) > mod1.fits Low Medium High Low 240.4943 189.1719 283.3337 High 326.5057 256.8281 384.6663 > tab1 Low Medium High Low 262 178 273 High 305 268 395 > mod1.fit[1,1]*mod1.fit[2,2]/(mod1.fit[1,2]*mod1.fit[2,1]) Error: object 'mod1.fit' not found > mod1.fit <- tapply(mod1$fitted, list(Cont,Sat), sum) > mod1.fit[1,1]*mod1.fit[2,2]/(mod1.fit[1,2]*mod1.fit[2,1]) [1] 1 > mosaicplot(mod1.fits, type="deviance", shade=T) > mod2 <- glm(Freq ~ Cont*Sat, family=poisson) > mod2 Call: glm(formula = Freq ~ Cont * Sat, family = poisson) Coefficients: (Intercept) ContHigh 3.08344 0.15197 SatMedium SatHigh -0.38656 0.04113 ContHigh:SatMedium ContHigh:SatHigh 0.25724 0.21745 Degrees of Freedom: 71 Total (i.e. Null); 66 Residual Null Deviance: 833.7 Residual Deviance: 745 AIC: 1098 > anova(mod1,mod2,test="LRT") Analysis of Deviance Table Model 1: Freq ~ Cont + Sat Model 2: Freq ~ Cont * Sat Resid. Df Resid. Dev Df Deviance Pr(>Chi) 1 68 750.17 2 66 745.04 2 5.1258 0.07708 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > chisq.test(tab1) Pearson's Chi-squared test data: tab1 X-squared = 5.1398, df = 2, p-value = 0.07654 > tab1 Low Medium High Low 262 178 273 High 305 268 395 > mod1.ll <- loglm(~1+2, data=tab1, param=TRUE, fit=TRUE) > mod1.ll Call: loglm(formula = ~1 + 2, data = tab1, param = TRUE, fit = TRUE) Statistics: X^2 df P(> X^2) Likelihood Ratio 5.125818 2 0.07708018 Pearson 5.139838 2 0.07654173 > mod1.ll$param $`(Intercept)` [1] 5.610201 $`1` Low High -0.1528753 0.1528753 $`2` Low Medium High 0.02537049 -0.21466986 0.18929936 > mod1$coef (Intercept) ContHigh SatMedium SatHigh 2.9977899 0.3057507 -0.2400404 0.1639289 > mod2.ll <- update(mod1.ll, .~.^2) > anova(mod1.ll,mod2.ll) LR tests for hierarchical log-linear models Model 1: ~1 + 2 Model 2: . ~ `1` + `2` + `1`:`2` Deviance df Delta(Dev) Delta(df) Model 1 5.125818 2 Model 2 0.000000 0 5.125818 2 Saturated 0.000000 0 0.000000 0 P(> Delta(Dev) Model 1 Model 2 0.07708 Saturated 1.00000 > tapply(mod1$fit,list(Cont,Sat),sum) # Same as 12*(linear predictor) Low Medium High Low 240.4943 189.1719 283.3337 High 326.5057 256.8281 384.6663 > mod1.ll$fit Low Medium High Low 240.4943 189.1719 283.3337 High 326.5057 256.8281 384.6663 > mod1 Call: glm(formula = Freq ~ Cont + Sat, family = poisson) Coefficients: (Intercept) ContHigh SatMedium SatHigh 2.9978 0.3058 -0.2400 0.1639 Degrees of Freedom: 71 Total (i.e. Null); 68 Residual Null Deviance: 833.7 Residual Deviance: 750.2 AIC: 1099 > options(contrasts=c("contr.sum","contr.poly")) # Sum to zero constraint > mod1b <- glm(Freq ~ Cont+Sat, family=poisson, data=housing) > mod1b Call: glm(formula = Freq ~ Cont + Sat, family = poisson, data = housing) Coefficients: (Intercept) Cont1 Sat1 Sat2 3.12529 -0.15288 0.02537 -0.21467 Degrees of Freedom: 71 Total (i.e. Null); 68 Residual Null Deviance: 833.7 Residual Deviance: 750.2 AIC: 1099 > options(contrasts=c("contr.treatment","contr.poly")) # Reference category equal to zero constraint > options()$contrasts [1] "contr.treatment" "contr.poly" >