Teng et al. (2020) analyzed the results of a survey of domestic cat owners in Australia. The survey focused on factors (e.g. cat demographics, owner attitudes and demographics, etc.) that might affect the prevalence of overweight and obese cats. They related nearly 1400 survey responses of owner-assessed body condition score [BCS with five categories: very underweight (1), somewhat underweight (2), ideal (3), chubby/overweight (4), and fat/obese (5)] to a range of categorical predictors with a multivariate multinomial GLM. We will use one aspect of their data to construct a contingency table relating the BCS, reduced to three categories (1&2, 3, 4&5) to cats’ begging behavior (four categories: never, sometimes, often, always).
The paper is here
Teng, K. T., McGreevy, P. D., Toribio, J. & Dhand, N. K. (2020). Positive attitudes towards feline obesity are strongly associated with ownership of obese cats. PLoS One, 15, e0234190.
Install the package epitools
library(tidyverse)
library(epitools)
library(vcd)Read in data
teng <- read_csv("../data/teng.csv")Rows: 12 Columns: 3── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (2): begging, bcs
dbl (1): count
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.tengMake begging a factor and put categories into a sensible order
teng$begging <- factor(teng$begging, levels=c("never","sometimes","often","always"))
teng.tab<-xtabs(count~begging + bcs, data=teng)
teng.tab bcs
begging bcs12 bcs3 bcs45
never 21 266 56
sometimes 55 527 162
often 14 106 69
always 10 47 44Generate chi-square statistic and look at residuals
chisq.test(teng.tab, correct=F)$exp bcs
begging bcs12 bcs3 bcs45
never 24.9092 235.641 82.450
sometimes 54.0305 511.129 178.841
often 13.7255 129.843 45.431
always 7.3348 69.387 24.278teng.chi<- chisq.test(teng.tab, correct=F)
teng.chi
Pearson's Chi-squared test
data: teng.tab
X-squared = 55.9, df = 6, p-value = 3e-10teng.chi$residuals bcs
begging bcs12 bcs3 bcs45
never -0.783268 1.977690 -2.912888
sometimes 0.131895 0.702024 -1.259312
often 0.074096 -2.092444 3.496684
always 0.984098 -2.687559 4.002581# Never vs sometimes begging
teng1 <- subset(teng,subset = begging %in% c('never','sometimes'))
teng1a <- subset(teng1,subset = bcs %in% c('bcs3','bcs45'))
teng1a.tab <- xtabs(count~begging + bcs, data=teng1a, drop.unused.levels=TRUE)
lodds(teng1a.tab, log = FALSE) odds for begging by bcs
bcs3 bcs45
0.50474 0.34568 loddsratio(teng1a.tab, log = FALSE) odds ratios for begging and bcs
[1] 1.4602# Never vs often
teng2 <- subset(teng,subset = begging %in% c('never','often'))
teng2a <- subset(teng2,subset = bcs %in% c('bcs3','bcs45'))
teng2a.tab <- xtabs(count~begging + bcs, data=teng2a, drop.unused.levels=TRUE)
lodds(teng2a.tab, log = FALSE) odds for begging by bcs
bcs3 bcs45
2.50943 0.81159 loddsratio(teng2a.tab, log = FALSE) odds ratios for begging and bcs
[1] 3.092# Never vs always begging
teng3 <- subset(teng,subset = begging %in% c('never','always'))
teng3a <- subset(teng3,subset = bcs %in% c('bcs3','bcs45'))
teng3a.tab <- xtabs(count~begging + bcs, data=teng3a, drop.unused.levels=TRUE)
lodds(teng3a.tab, log = FALSE) odds for begging by bcs
bcs3 bcs45
5.6596 1.2727 loddsratio(teng3a.tab, log = FALSE) odds ratios for begging and bcs
[1] 4.4468tengfull.glm <- glm(count~begging*bcs, family=poisson, data=teng)
summary(tengfull.glm)
Call:
glm(formula = count ~ begging * bcs, family = poisson, data = teng)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.1241 0.0482 85.52 < 2e-16 ***
begging1 0.0937 0.0786 1.19 0.2330
begging2 0.9966 0.0615 16.20 < 2e-16 ***
begging3 -0.2786 0.0873 -3.19 0.0014 **
bcs1 -1.1257 0.0853 -13.20 < 2e-16 ***
bcs2 0.9670 0.0555 17.42 < 2e-16 ***
begging1:bcs1 -0.0476 0.1380 -0.34 0.7304
begging2:bcs1 0.0123 0.1084 0.11 0.9094
begging3:bcs1 -0.0808 0.1564 -0.52 0.6056
begging1:bcs2 0.3987 0.0869 4.59 4.5e-06 ***
begging2:bcs2 0.1795 0.0697 2.58 0.0100 *
begging3:bcs2 -0.1491 0.0998 -1.49 0.1351
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 1.5702e+03 on 11 degrees of freedom
Residual deviance: 3.2419e-14 on 0 degrees of freedom
AIC: 95.6
Number of Fisher Scoring iterations: 3anova(tengfull.glm, test="Chisq")Analysis of Deviance Table
Model: poisson, link: log
Response: count
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL 11 1570
begging 3 670 8 900 < 2e-16 ***
bcs 2 847 6 53 < 2e-16 ***
begging:bcs 6 53 0 0 1.1e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1tengred.glm <- glm(count~begging+bcs, family=poisson, data=teng)
summary(tengred.glm)
Call:
glm(formula = count ~ begging + bcs, family = poisson, data = teng)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.1022 0.0453 90.52 < 2e-16 ***
begging1 0.2611 0.0517 5.05 4.4e-07 ***
begging2 1.0354 0.0434 23.83 < 2e-16 ***
begging3 -0.3349 0.0621 -5.39 7.1e-08 ***
bcs1 -1.1480 0.0700 -16.40 < 2e-16 ***
bcs2 1.0991 0.0438 25.10 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 1570.221 on 11 degrees of freedom
Residual deviance: 53.239 on 6 degrees of freedom
AIC: 136.8
Number of Fisher Scoring iterations: 4anova(tengred.glm, test="Chisq")Analysis of Deviance Table
Model: poisson, link: log
Response: count
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL 11 1570
begging 3 670 8 900 <2e-16 ***
bcs 2 847 6 53 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1anova(tengfull.glm, tengred.glm, test="Chisq")Analysis of Deviance Table
Model 1: count ~ begging * bcs
Model 2: count ~ begging + bcs
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 0 0.0
2 6 53.2 -6 -53.2 1.1e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1lrtest(tengfull.glm, tengred.glm)Likelihood ratio test
Model 1: count ~ begging * bcs
Model 2: count ~ begging + bcs
#Df LogLik Df Chisq Pr(>Chisq)
1 12 -35.8
2 6 -62.4 -6 53.2 1.1e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1