################################ #### Example: Vitamin C & French skiers #### LOGLIN() & GLM() functions #### Related: See VitaminC.R ################################# #### Here is one way to read the data vector of values with labels for the table ski<-matrix(c(31, 17, 109, 122), ncol=2, dimnames=list(Treatment=c("Placebo", "VitaminC"), Cold=c("Cold", "NoCold"))) ski #### Here is how we did this previously, for more details see VitaminC.R #### Pearson's Chi-squared test with Yates' continuity correction result<-chisq.test(ski) result ###################### #### VIA LOGLIN() #### #### Fitting 2-way log-linear model of independence #### using the table structure and loglin() function; this is similar to CATMOD in SAS #### first check if the data are in the table format is.table(ski) #### if not make them to be in table format with ski<-as.table(ski) #### fit the model #### list(1,2) indicates that we are fitting the two margins independently #### fit=TRUE and param=TRUE, say that we want have the fitted values of the cells #### and the estimated model parameters returned #### for other options see R help, e.g., type the ?loglin in the prompt ski.ind<-loglin(ski, list(1, 2), fit=TRUE, param=TRUE) ski.ind #### get the p-value; note that this output is the same as the chi-squared test #### WITHOUT Yates' continuity correction. 1-pchisq(ski.ind\$lrt, ski.ind\$df) #### Fitting 2-way saturated log-linear model #### using the table structure and loglin() function; ####this is similar to CATMOD in SAS ski.sat<-loglin(ski, list(c(1, 2)), fit=TRUE, param=TRUE) #### notice the zero value for LRT, and the PERFECT fit!!! ski.sat 1-pchisq(ski.sat\$lrt, ski.sat\$df) ################### #### VIA GLM() #### #### Fitting 2-way log-linear model of independence #### using the database structure and glm() function; this is similar to GENMOD in SAS #### first we must format the data into a data frame format ski.data<-as.data.frame(ski) ski.data #### using the glm() to fit the loglinear model is similar to fitting a regression model #### but the trick is in specifying what the correct response, and the family of distributions #### (i.e., sampling scheme!) #### we need to specify the response that we are modeling which are the cell counts; #### ski.data\$Freq #### this line: ski.data\$Treatment+ski.data\$Cold, specifies that we have two main effects, #### and no interaction term #### we need to specify that sampling distribution, that is for log-linear models when modeling #### the counts we assume the most general random mechanism of Poisson sampling: family=poinsson(). #### This assures that errors follow Poisson model AND that we are modeling "log" of the response #### that is of our counts. ski.ind<-glm(ski.data\$Freq~ski.data\$Treatment+ski.data\$Cold, family=poisson()) #### to view the model and the relevant statistics ski.ind #### another way of viewing the model summary(ski.ind) #### this is ANOVA table similar to what you have seen in linear regression models that looks #### at the contribution of different model parameters to the overall fit of the model anova(ski.ind) #### looking at and saving the fitted values of the cell counts fits<-fitted(ski.ind) #### looking at and saving the pearson residuals resids <- residuals(ski.ind,type="pearson") #### checking for the influential points h <- lm.influence(ski.ind)\$hat #### adjusted pearson residuals, e.g., standardized adjresids <- resids/sqrt(1-h) #### putting them together into a table round(cbind(ski.data\$Freq,fits,adjresids),2) #### Fitting 2-way saturated log-linear model #### using the database structure and glm() function; this is similar to GENMOD in SAS #### Notice now that we have an interaction term in the model: ski.data\$Treatment*ski.data\$Cold ski.sat<-glm(ski.data\$Freq~ski.data\$Treatment*ski.data\$Cold, family=poisson()) ski.sat anova(ski.sat)