code 9.1

vote<-c("yes","yes","yes","no","no","no","undec","undec","undec")
urban<-c("urban","sub","rural","urban","sub","rural","urban","sub","rural")
count<-c(8, 17, 7, 6, 8, 15, 19, 7, 11)

vote.data<-data.frame(vote, urban,  count)

code 9.2

vote.table<-xtabs(count ~ vote + urban)

vote.table

##        urban
## vote    rural sub urban
##   no       15   8     6
##   undec    11   7    19
##   yes       7  17     8

apply(vote.table, 1, sum)  #행합

##    no undec   yes 
##    29    37    32

apply(vote.table, 2, sum)   #열합

## rural   sub urban 
##    33    32    33

barplot(vote.table)     #누적막대그림

barplot(vote.table, beside=TRUE, legend=TRUE, ylim=c(0,20)) #막대그림

code 9.3

install.packages(“Deducer”) install.packages(“vcd”) 이것들 필요함

# code 9.2의 자료를 사용

chisq.test(vote.table)$expected #기대빈도의 생성

##        urban
## vote        rural       sub     urban
##   no     9.765306  9.469388  9.765306
##   undec 12.459184 12.081633 12.459184
##   yes   10.775510 10.448980 10.775510

summary(vote.table)     #카이제곱 검정

## Call: xtabs(formula = count ~ vote + urban)
## Number of cases in table: 98 
## Number of factors: 2 
## Test for independence of all factors:
##  Chisq = 16.373, df = 4, p-value = 0.002558

chisq.test(vote.table)  #카이제곱 검정

## 
##  Pearson's Chi-squared test
## 
## data:  vote.table
## X-squared = 16.373, df = 4, p-value = 0.002558

library(Deducer)

## Loading required package: ggplot2
## Loading required package: JGR
## Loading required package: rJava
## Loading required package: JavaGD
## Loading required package: iplots
## 
## Please type JGR() to launch console. Platform specific launchers (.exe and .app) can also be obtained at http://www.rforge.net/JGR/files/.
## 
## 
## Loading required package: car
## Loading required package: MASS
## 
## 
## Note Non-JGR console detected:
##  Deducer is best used from within JGR (http://jgr.markushelbig.org/).
##  To Bring up GUI dialogs, type deducer().
## 
## 
## Attaching package: 'Deducer'
## 
## The following object is masked from 'package:stats':
## 
##     summary.lm

likelihood.test(vote.table)     #우도비 카이제곱

## 
##  Log likelihood ratio (G-test) test of independence without
##  correction
## 
## data:  vote.table
## Log likelihood ratio statistic (G) = 15.731, X-squared df = 4,
## p-value = 0.003402

library(vcd)

## Loading required package: grid

assocstats(vote.table)      #연관성 분석을 통해 얻어진 통계량

##                     X^2 df  P(> X^2)
## Likelihood Ratio 15.731  4 0.0034019
## Pearson          16.373  4 0.0025575
## 
## Phi-Coefficient   : NA 
## Contingency Coeff.: 0.378 
## Cramer's V        : 0.289

vote.data2<-c()

for(i in 1:nrow(vote.data))
{
    for(j in 1:vote.data[i,3])
    {
        vote.data2<-rbind(vote.data2, vote.data[i,1:2])
    }
}

head(vote.data2)

##   vote urban
## 1  yes urban
## 2  yes urban
## 3  yes urban
## 4  yes urban
## 5  yes urban
## 6  yes urban

ordered(c("no","undec", "yes"))     #서열형 변수 설정

## [1] no    undec yes  
## Levels: no < undec < yes

head(vote.data2)

##   vote urban
## 1  yes urban
## 2  yes urban
## 3  yes urban
## 4  yes urban
## 5  yes urban
## 6  yes urban

table(vote.data2$vote, vote.data2$urban)

##        
##         rural sub urban
##   no       15   8     6
##   undec    11   7    19
##   yes       7  17     8

vote.data2$vote<-factor(vote.data2$vote, ordered=T)
vote.data2$vote

##  [1] yes   yes   yes   yes   yes   yes   yes   yes   yes   yes   yes  
## [12] yes   yes   yes   yes   yes   yes   yes   yes   yes   yes   yes  
## [23] yes   yes   yes   yes   yes   yes   yes   yes   yes   yes   no   
## [34] no    no    no    no    no    no    no    no    no    no    no   
## [45] no    no    no    no    no    no    no    no    no    no    no   
## [56] no    no    no    no    no    no    undec undec undec undec undec
## [67] undec undec undec undec undec undec undec undec undec undec undec
## [78] undec undec undec undec undec undec undec undec undec undec undec
## [89] undec undec undec undec undec undec undec undec undec undec
## Levels: no < undec < yes

vote.data2$urban<-factor(vote.data2$urban, ordered=T)

ordered(c("rural", "sub", "urban")) #서열형 변수 설정

## [1] rural sub   urban
## Levels: rural < sub < urban

v1<-as.numeric(vote.data2$vote) #서열형 변수를 숫자로 변환
u1<-as.numeric(vote.data2$urban)    #서열형 변수를 숫자로 변환

data.frame(v1, u1)

##    v1 u1
## 1   3  3
## 2   3  3
## 3   3  3
## 4   3  3
## 5   3  3
## 6   3  3
## 7   3  3
## 8   3  3
## 9   3  2
## 10  3  2
## 11  3  2
## 12  3  2
## 13  3  2
## 14  3  2
## 15  3  2
## 16  3  2
## 17  3  2
## 18  3  2
## 19  3  2
## 20  3  2
## 21  3  2
## 22  3  2
## 23  3  2
## 24  3  2
## 25  3  2
## 26  3  1
## 27  3  1
## 28  3  1
## 29  3  1
## 30  3  1
## 31  3  1
## 32  3  1
## 33  1  3
## 34  1  3
## 35  1  3
## 36  1  3
## 37  1  3
## 38  1  3
## 39  1  2
## 40  1  2
## 41  1  2
## 42  1  2
## 43  1  2
## 44  1  2
## 45  1  2
## 46  1  2
## 47  1  1
## 48  1  1
## 49  1  1
## 50  1  1
## 51  1  1
## 52  1  1
## 53  1  1
## 54  1  1
## 55  1  1
## 56  1  1
## 57  1  1
## 58  1  1
## 59  1  1
## 60  1  1
## 61  1  1
## 62  2  3
## 63  2  3
## 64  2  3
## 65  2  3
## 66  2  3
## 67  2  3
## 68  2  3
## 69  2  3
## 70  2  3
## 71  2  3
## 72  2  3
## 73  2  3
## 74  2  3
## 75  2  3
## 76  2  3
## 77  2  3
## 78  2  3
## 79  2  3
## 80  2  3
## 81  2  2
## 82  2  2
## 83  2  2
## 84  2  2
## 85  2  2
## 86  2  2
## 87  2  2
## 88  2  1
## 89  2  1
## 90  2  1
## 91  2  1
## 92  2  1
## 93  2  1
## 94  2  1
## 95  2  1
## 96  2  1
## 97  2  1
## 98  2  1

?cor.test

## starting httpd help server ... done

cor.test(v1, u1, method="kendall", exact=FALSE)     #켄달 타우-b

## 
##  Kendall's rank correlation tau
## 
## data:  v1 and u1
## z = 1.4251, p-value = 0.1541
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##       tau 
## 0.1287195

code 9. 4

p.value<-1-pchisq(16.3729, df=4)

chi<-seq(10, 20, by=0.01)
pdf<-dchisq(chi, df=4)

chi.graph<-data.frame(chi, pdf)
chi.graph$pdf2<-with(chi.graph, ifelse(chi>=16.3729, pdf, NA))

chi.area<-data.frame(chi.graph$chi, chi.graph$pdf2)
names(chi.area) <- c("chi", "pdf")
head(chi.area)

##     chi pdf
## 1 10.00  NA
## 2 10.01  NA
## 3 10.02  NA
## 4 10.03  NA
## 5 10.04  NA
## 6 10.05  NA

ggplot(chi.graph, aes(x=chi, y=pdf)) + geom_line() +
    geom_area(data=chi.area, alpha=0.3)

## Warning: Removed 638 rows containing missing values (position_stack).

code 9.5

welfare<-c("plural", "plural", "flexible", "flexible")
gender<-c("m", "f", "m", "f")
count<-c(12, 3, 3, 14)

small<-data.frame(welfare, gender, count)

small.table<-xtabs(count ~ welfare + gender)    #정확성검정을 실시할 테이블 작성
small.table 

##           gender
## welfare     f  m
##   flexible 14  3
##   plural    3 12

fisher.test(small.table)    #피셔의 정확성검정

## 
##  Fisher's Exact Test for Count Data
## 
## data:  small.table
## p-value = 0.001033
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##    2.507119 159.407427
## sample estimates:
## odds ratio 
##   16.41917

code 9.6

region<-c("urban","urban","urban","urban","rural","rural","rural","rural")
inequality<-c("low","low","high","high","low","low","high","high")
trust<-c("good","bad","good","bad","good","bad","good","bad")
count<-c(48, 12, 96, 94, 55, 130, 7, 53)

inequ<-data.frame(region, inequality, trust, count)

reg.false<-xtabs(count ~ inequality + trust, data=inequ)

#list형식으로 데이터를 만든다. 이래야 CMH통계량 산출이 가능하다
reg.true<-xtabs(count ~ inequality + trust + region , data=inequ)

chisq.test(reg.false, correct=F)    ##2x2 table에서 correction

## 
##  Pearson's Chi-squared test
## 
## data:  reg.false
## X-squared = 0.036004, df = 1, p-value = 0.8495

chisq.test(reg.false)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  reg.false
## X-squared = 0.0097136, df = 1, p-value = 0.9215

summary(reg.false)

## Call: xtabs(formula = count ~ inequality + trust, data = inequ)
## Number of cases in table: 495 
## Number of factors: 2 
## Test for independence of all factors:
##  Chisq = 0.036, df = 1, p-value = 0.8495

#CMH통계량의 산출
mantelhaen.test(reg.true, correct=FALSE)    # 연속성 수정

## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  reg.true
## Mantel-Haenszel X-squared = 23.639, df = 1, p-value = 1.162e-06
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  2.099716 6.166385
## sample estimates:
## common odds ratio 
##          3.598285

mantelhaen.test(reg.true)

## 
##  Mantel-Haenszel chi-squared test with continuity correction
## 
## data:  reg.true
## Mantel-Haenszel X-squared = 22.558, df = 1, p-value = 2.039e-06
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  2.099716 6.166385
## sample estimates:
## common odds ratio 
##          3.598285

code 9.7

before<-c("yes", "yes", "no", "no")
after<-c("yes", "no", "yes", "no")
count<-c(19, 5, 15, 48)

road<-data.frame(before, after, count)

road.table<-xtabs(count ~ before + after)

#맥니마 테스트
mcnemar.test(road.table, correct = FALSE)   # 연속성 수정(?)

## 
##  McNemar's Chi-squared test
## 
## data:  road.table
## McNemar's chi-squared = 5, df = 1, p-value = 0.02535

mcnemar.test(road.table)

## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  road.table
## McNemar's chi-squared = 4.05, df = 1, p-value = 0.04417

code 9.8

#CMH통계량을 산출하기 위해서 list형식으로 데이터를 수정해야 한다. 
road

##   before after count
## 1    yes   yes    19
## 2    yes    no     5
## 3     no   yes    15
## 4     no    no    48

road2<-c()

temp1<-c()
temp2<-c()
temp<-c()


for(j in 1:2)
{
    id<-0
    temp3<-
    for(i in 1:4)
    {
        for(k in 1:road[i,3])
        {
            temp1<-road[i,j]
            temp2<-names(road)[j]
            id<-id+1
            temp<-data.frame(id, temp1, temp2)
            road2<-rbind(road2, temp)
        }
    }
}

road2

##     id temp1  temp2
## 1    1   yes before
## 2    2   yes before
## 3    3   yes before
## 4    4   yes before
## 5    5   yes before
## 6    6   yes before
## 7    7   yes before
## 8    8   yes before
## 9    9   yes before
## 10  10   yes before
## 11  11   yes before
## 12  12   yes before
## 13  13   yes before
## 14  14   yes before
## 15  15   yes before
## 16  16   yes before
## 17  17   yes before
## 18  18   yes before
## 19  19   yes before
## 20  20   yes before
## 21  21   yes before
## 22  22   yes before
## 23  23   yes before
## 24  24   yes before
## 25  25    no before
## 26  26    no before
## 27  27    no before
## 28  28    no before
## 29  29    no before
## 30  30    no before
## 31  31    no before
## 32  32    no before
## 33  33    no before
## 34  34    no before
## 35  35    no before
## 36  36    no before
## 37  37    no before
## 38  38    no before
## 39  39    no before
## 40  40    no before
## 41  41    no before
## 42  42    no before
## 43  43    no before
## 44  44    no before
## 45  45    no before
## 46  46    no before
## 47  47    no before
## 48  48    no before
## 49  49    no before
## 50  50    no before
## 51  51    no before
## 52  52    no before
## 53  53    no before
## 54  54    no before
## 55  55    no before
## 56  56    no before
## 57  57    no before
## 58  58    no before
## 59  59    no before
## 60  60    no before
## 61  61    no before
## 62  62    no before
## 63  63    no before
## 64  64    no before
## 65  65    no before
## 66  66    no before
## 67  67    no before
## 68  68    no before
## 69  69    no before
## 70  70    no before
## 71  71    no before
## 72  72    no before
## 73  73    no before
## 74  74    no before
## 75  75    no before
## 76  76    no before
## 77  77    no before
## 78  78    no before
## 79  79    no before
## 80  80    no before
## 81  81    no before
## 82  82    no before
## 83  83    no before
## 84  84    no before
## 85  85    no before
## 86  86    no before
## 87  87    no before
## 88   1   yes  after
## 89   2   yes  after
## 90   3   yes  after
## 91   4   yes  after
## 92   5   yes  after
## 93   6   yes  after
## 94   7   yes  after
## 95   8   yes  after
## 96   9   yes  after
## 97  10   yes  after
## 98  11   yes  after
## 99  12   yes  after
## 100 13   yes  after
## 101 14   yes  after
## 102 15   yes  after
## 103 16   yes  after
## 104 17   yes  after
## 105 18   yes  after
## 106 19   yes  after
## 107 20    no  after
## 108 21    no  after
## 109 22    no  after
## 110 23    no  after
## 111 24    no  after
## 112 25   yes  after
## 113 26   yes  after
## 114 27   yes  after
## 115 28   yes  after
## 116 29   yes  after
## 117 30   yes  after
## 118 31   yes  after
## 119 32   yes  after
## 120 33   yes  after
## 121 34   yes  after
## 122 35   yes  after
## 123 36   yes  after
## 124 37   yes  after
## 125 38   yes  after
## 126 39   yes  after
## 127 40    no  after
## 128 41    no  after
## 129 42    no  after
## 130 43    no  after
## 131 44    no  after
## 132 45    no  after
## 133 46    no  after
## 134 47    no  after
## 135 48    no  after
## 136 49    no  after
## 137 50    no  after
## 138 51    no  after
## 139 52    no  after
## 140 53    no  after
## 141 54    no  after
## 142 55    no  after
## 143 56    no  after
## 144 57    no  after
## 145 58    no  after
## 146 59    no  after
## 147 60    no  after
## 148 61    no  after
## 149 62    no  after
## 150 63    no  after
## 151 64    no  after
## 152 65    no  after
## 153 66    no  after
## 154 67    no  after
## 155 68    no  after
## 156 69    no  after
## 157 70    no  after
## 158 71    no  after
## 159 72    no  after
## 160 73    no  after
## 161 74    no  after
## 162 75    no  after
## 163 76    no  after
## 164 77    no  after
## 165 78    no  after
## 166 79    no  after
## 167 80    no  after
## 168 81    no  after
## 169 82    no  after
## 170 83    no  after
## 171 84    no  after
## 172 85    no  after
## 173 86    no  after
## 174 87    no  after

road2$count<-1      #테이블을 구성하기 위해 카운트 정보 변수 생성
head(road2)

##   id temp1  temp2 count
## 1  1   yes before     1
## 2  2   yes before     1
## 3  3   yes before     1
## 4  4   yes before     1
## 5  5   yes before     1
## 6  6   yes before     1

road2.table<-xtabs(count ~ temp1 + temp2 + id, data=road2)

#CMH통계량 산출. 맥니마 테스트의 결과와 같다.
mantelhaen.test(road2.table, correct=FALSE)

## 
##  Mantel-Haenszel chi-squared test without continuity correction
## 
## data:  road2.table
## Mantel-Haenszel X-squared = 5, df = 1, p-value = 0.02535
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.090342 8.254292
## sample estimates:
## common odds ratio 
##                 3

mantelhaen.test(road2.table)

## 
##  Mantel-Haenszel chi-squared test with continuity correction
## 
## data:  road2.table
## Mantel-Haenszel X-squared = 4.05, df = 1, p-value = 0.04417
## alternative hypothesis: true common odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.090342 8.254292
## sample estimates:
## common odds ratio 
##                 3