library(vcd)
## Warning: package 'vcd' was built under R version 3.5.1
## Loading required package: grid
library(vcdExtra)
## Warning: package 'vcdExtra' was built under R version 3.5.1
## Loading required package: gnm
library(ca)
## Warning: package 'ca' was built under R version 3.5.1
library(logmult)
## Warning: package 'logmult' was built under R version 3.5.1
##
## Attaching package: 'logmult'
## The following object is masked from 'package:gnm':
##
## se
## The following object is masked from 'package:vcd':
##
## assoc
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.5.1
Exercise 6.2
data("criminal", package = "logmult")
criminal
## Age
## Year 15 16 17 18 19
## 1955 141 285 320 441 427
## 1956 144 292 342 441 396
## 1957 196 380 424 462 427
## 1958 212 424 399 442 430
(a) What percentages of Pearson chisq for association are explained by the various dimensions?
crim <- margin.table(criminal, 1:2)
(crim.ca <- ca(crim))
##
## Principal inertias (eigenvalues):
## 1 2 3
## Value 0.004939 0.000491 3.8e-05
## Percentage 90.33% 8.98% 0.69%
##
##
## Rows:
## 1955 1956 1957 1958
## Mass 0.229751 0.229893 0.268897 0.271459
## ChiDist 0.090897 0.061048 0.047585 0.088033
## Inertia 0.001898 0.000857 0.000609 0.002104
## Dim. 1 1.253085 0.827543 -0.553684 -1.212927
## Dim. 2 -0.984738 0.733468 1.206411 -0.982745
##
##
## Columns:
## 15 16 17 18 19
## Mass 0.098648 0.196584 0.211388 0.254235 0.239146
## ChiDist 0.101134 0.093089 0.044072 0.071068 0.066594
## Inertia 0.001009 0.001703 0.000411 0.001284 0.001061
## Dim. 1 -1.433374 -1.297270 -0.332608 1.000960 0.887539
## Dim. 2 -0.333181 -0.808352 1.676250 0.307874 -1.007063
(a) ans: Dim 1(90.3%) and Dim 2(8.98%)
(b) Plot the 2D correspondence analysis solution. Describe the pattern of association between year and age.
plot(crim.ca)
(b) ans: (1) There seems to be a slight association between age 17 and the year 1957. (2) There seems to be a strong association between age 16 and the year 1958. (3) There is a slight association between age 18 and the year 1956.
Exercise 6.11
data("Vietnam", package = "vcdExtra")
str(Vietnam)
## 'data.frame': 40 obs. of 4 variables:
## $ sex : Factor w/ 2 levels "Female","Male": 1 1 1 1 1 1 1 1 1 1 ...
## $ year : int 1 1 1 1 2 2 2 2 3 3 ...
## $ response: Factor w/ 4 levels "A","B","C","D": 1 2 3 4 1 2 3 4 1 2 ...
## $ Freq : int 13 19 40 5 5 9 33 3 22 29 ...
(a)
Vietnam$yearsex<- paste(Vietnam$year, Vietnam$sex, sep = ':')
VN <- xtabs(Freq ~ response + Vietnam$yearsex, data = Vietnam)
viet <- margin.table(VN, 1:2)
viet.ca <- ca(viet)
summary(viet.ca)
##
## Principal inertias (eigenvalues):
##
## dim value % cum% scree plot
## 1 0.085680 73.6 73.6 ******************
## 2 0.027881 23.9 97.5 ******
## 3 0.002854 2.5 100.0 *
## -------- -----
## Total: 0.116415 100.0
##
##
## Rows:
## name mass qlt inr k=1 cor ctr k=2 cor ctr
## 1 | A | 255 985 381 | -414 985 509 | 1 0 0 |
## 2 | B | 235 720 60 | -135 608 50 | -58 112 28 |
## 3 | C | 419 999 283 | 247 773 298 | 133 226 267 |
## 4 | D | 92 995 276 | 366 383 143 | -463 612 705 |
##
## Columns:
## name mass qlt inr k=1 cor ctr k=2 cor ctr
## 1 | 1Fml | 24 818 13 | 167 452 8 | 150 367 20 |
## 2 | 1Mal | 139 997 181 | -386 986 242 | 41 11 8 |
## 3 | 2Fml | 16 995 35 | 407 647 31 | 299 349 51 |
## 4 | 2Mal | 140 984 131 | -326 982 175 | 15 2 1 |
## 5 | 3Fml | 53 999 112 | 334 453 69 | 367 547 256 |
## 6 | 3Mal | 138 904 40 | -175 904 49 | 4 0 0 |
## 7 | 4Fml | 32 982 37 | 344 887 44 | 113 95 15 |
## 8 | 4Mal | 149 383 23 | -81 372 11 | -14 11 1 |
## 9 | 5Fml | 59 994 153 | 453 686 143 | 304 309 197 |
## 10 | 5Mal | 248 1000 276 | 281 608 228 | -225 391 451 |
A = Defeat Vietnam by widespread Bombing:
There is a somewhat strong association with both male and female here. Greater for male than for female.
B = Maintain present policy:
There is a somewhat strong association with both male and female here.
C = De-escalate military activity: stop bombing, … :
There is an association with both male and female.
D = Withdraw military forces:
No association. A weak association with males.
(c) Use mjca() to carry out an MCA on the three-way table. Make a useful plot of the solution and interpret in terms of the relationship of the response to year and sex.
data("Vietnam", package = "vcdExtra")
viet2 <- xtabs(Freq ~ sex + year + response, data = Vietnam)
data(viet2)
## Warning in data(viet2): data set 'viet2' not found
viet.mca <- mjca(viet2)
summary(viet.mca)
##
## Principal inertias (eigenvalues):
##
## dim value % cum% scree plot
## 1 0.028219 65.1 65.1 ********************
## 2 0.007445 17.2 82.3 *****
## 3 0.000380 0.9 83.2
## 4 1e-06000 0.0 83.2
## -------- -----
## Total: 0.043317
##
##
## Columns:
## name mass qlt inr k=1 cor ctr k=2 cor ctr
## 1 | sex:Female | 62 853 105 | 328 749 235 | -122 104 123 |
## 2 | sex:Male | 272 853 24 | -74 749 53 | 28 104 28 |
## 3 | year:1 | 55 856 103 | -201 780 78 | -63 76 29 |
## 4 | year:2 | 52 948 104 | -231 948 99 | -6 1 0 |
## 5 | year:3 | 64 899 100 | 104 306 24 | -145 592 179 |
## 6 | year:4 | 60 55 99 | -12 44 0 | 6 10 0 |
## 7 | year:5 | 102 756 88 | 167 490 102 | 123 267 210 |
## 8 | response:A | 85 818 95 | -260 797 203 | -42 21 20 |
## 9 | response:B | 78 663 93 | -93 662 24 | -4 1 0 |
## 10 | response:C | 140 856 74 | 182 817 164 | -39 38 29 |
## 11 | response:D | 31 736 113 | 125 107 17 | 304 629 381 |
plot(viet.mca)
Interpretation:
A = Defeat Vietnam by widespread Bombing:
There is an association between Year 1 and Year 2
B = Maintain present policy:
There is an association between Males and Year 4
C = De-escalate military activity: stop bombing, … :
There is virtually no association in terms of the relationship
of Response C to year and sex
D = Withdraw military forces:
There is no association in terms of the relationship of Response D to year and sex.
