Exercise 4.2 The data set Abortion in vcdExtra gives a 2 ? 2 ? 2 table of opinions regarding abortion in relation to sex and status of the respondent.
data("Abortion", package = "vcdExtra")
str(Abortion)
## table [1:2, 1:2, 1:2] 171 152 138 167 79 148 112 133
## - attr(*, "dimnames")=List of 3
## ..$ Sex : chr [1:2] "Female" "Male"
## ..$ Status : chr [1:2] "Lo" "Hi"
## ..$ Support_Abortion: chr [1:2] "Yes" "No"
## Loading required package: grid
## Status Lo Hi
## Sex Support_Abortion
## Female Yes 171 138
## No 79 112
## Male Yes 152 167
## No 148 133
Abortion3<-aperm(Abortion, c(2,3,1))
fourfold(Abortion3)
# Sex by support for abortion, stratified by status
summary(oddsratio(Abortion2))
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## Lo 0.74555 0.17844 4.1781 2.94e-05 ***
## Hi -0.01889 0.17228 -0.1096 0.9127
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Status by support for abortion, stratified by sex
summary(oddsratio(Abortion3))
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## Female 0.56346 0.18623 3.0256 0.002481 **
## Male -0.20098 0.16384 -1.2267 0.219941
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Exercise 4.4 The Hospital data in vcd gives a 3 ? 3 table relating the length of stay (in years) of 132 long-term schizophrenic patients in two London mental hospitals with the frequency of visits by family and friends.
(a)Carry out a ??2 test for association between the two variables.
data("Hospital", package="vcd")
chisq.test(Hospital)
##
## Pearson's Chi-squared test
##
## data: Hospital
## X-squared = 35.171, df = 4, p-value = 4.284e-07
(b)Use assocstats () to compute association statistics. How would you describe the strength of association here?
assocstats(Hospital)
## X^2 df P(> X^2)
## Likelihood Ratio 38.353 4 9.4755e-08
## Pearson 35.171 4 4.2842e-07
##
## Phi-Coefficient : NA
## Contingency Coeff.: 0.459
## Cramer's V : 0.365
(c)Produce an association plot for these data, with visit frequency as the vertical variable. Describe the pattern of the relation you see here.
assoc(Hospital, shade=TRUE)
(d)Both variables can be considered ordinal, so CMHtest () may be useful here. Carry out that analysis. Do any of the tests lead to different conclusions?
library(vcdExtra)
## Loading required package: gnm
## Warning: package 'gnm' was built under R version 3.4.4
CMHtest(Hospital)
## Cochran-Mantel-Haenszel Statistics for Visit frequency by Length of stay
##
## AltHypothesis Chisq Df Prob
## cor Nonzero correlation 29.138 1 6.7393e-08
## rmeans Row mean scores differ 34.391 2 3.4044e-08
## cmeans Col mean scores differ 29.607 2 3.7233e-07
## general General association 34.905 4 4.8596e-07
Exercise 4.6 The two-way table Mammograms in vcdExtra gives ratings on the severity of diagnosis of 110 mammograms by two raters.
(a)Assess the strength of agreement between the raters using Cohen’s ??, both unweighted and weighted.
Kappa(Mammograms)
## value ASE z Pr(>|z|)
## Unweighted 0.3713 0.06033 6.154 7.560e-10
## Weighted 0.5964 0.04923 12.114 8.901e-34
Kappa(Mammograms, weights= "Fleiss-Cohen")
## value ASE z Pr(>|z|)
## Unweighted 0.3713 0.06033 6.154 7.560e-10
## Weighted 0.7641 0.03996 19.122 1.667e-81
(b)Use agreementplot () for a graphical display of agreement here.
agreementplot(Mammograms, main="Unweighted", weights=1)
agreementplot(Mammograms, main="Weighted")
(c)Compare the Kappa measures with the results from assocstats (). What is a reasonable interpretation of each of these measures
assocstats(Mammograms)
## X^2 df P(> X^2)
## Likelihood Ratio 92.619 9 4.4409e-16
## Pearson 83.516 9 3.2307e-14
##
## Phi-Coefficient : NA
## Contingency Coeff.: 0.657
## Cramer's V : 0.503
Exercise 4.7 Agresti and Winner (1997) gave the data in below on the ratings of 160 movies by the reviewers Gene Siskel and Roger Ebert for the period from April 1995 through September 1996. The rating categories were Con (“thumbs down”), Mixed, and Pro (“thumbs up”).
ratings <- matrix(
c(
24, 8, 13,
8, 13, 11,
10, 9, 64 ), 3, 3, byrow=TRUE)
dimnames(ratings) <- list(Siskal=c("Con", "Mixed", "Pro"),
Ebert =c("Con", "Mixed", "Pro"))
Kappa(ratings)
## value ASE z Pr(>|z|)
## Unweighted 0.3888 0.05979 6.503 7.870e-11
## Weighted 0.4269 0.06350 6.723 1.781e-11
agreementplot(ratings)