R Correlation Coefficient of Categorical Variables (カテゴリカルデータの相関係数)

Program R

#全ての説明変数が量的であれば、主成分分析、因子分析等の

#量的変数を対象とした手法の適用が可能。

#カテゴリカル変数が説明変数に含まれている場合に、

#量的変数を対象とした手法の適用が出来るようにする方法を

#下記で確認する。

 

#In case of ordinal categorical variables

#カテゴリカル説明変数が、本来連続的な変数であると考え、

#その連続的な変数を閾値によって、カテゴリカルに分断している

#と考える。

 

#連続変数の分布を仮定し、最尤推定によってパラメータを決定。

#決まった相関係数をテトラコリック相関(カテゴリ数が2の場合)

#／ポリコリック相関(カテゴリ数が3以上の場合)という

#Rではpolycor,hectorを用いることで計算可能

 

library(vcd)

head(DanishWelfare)

 

tabledata <- xtabs(Freq~Alcohol+Income,data=DanishWelfare)

tabledata

 

install.packages("polycor")

library(polycor)

 

polychor(tabledata,ML=TRUE)

#MLは最尤推定法を示す。Maximum-Likelihood

#付さないと2段階推定法という簡便な方法。

 

#2変数はテーブルではなくペアで与えられていても良い。

x <- c(1,1,2,2,1,3,3,3,1,1)

y <- c(3,3,1,1,2,2,2,1,2,2)

polychor(x,y,ML=TRUE)

 

#個票データでも良い。集計データの場合は一度個票かテーブルに

#変換する必要がある。

 

#仕組みを確認

#標準正規分布に従う100個の乱数列r1,r2を準備、変数z,wを作成

set.seed(0)

r1 <- rnorm(100)

r2 <- rnorm(100)

z <- r1+r2

w <- r1

cor(z, w)

 

#仮定している正規分布との対応を確認

x <- cut(z, breaks <- c(-100,-2,-1,0,1,2,100), label<-c(1:6),

         right = FALSE)

y <- cut(w, breaks <- c(-100,-2,-1,0,1,2,100), label<-c(1:6),

         right = FALSE)

polychor(x,y)

 

#########順序カテゴリカル変数と連続変数との相関#########

library(polycor)

x <- c(1,7,2,5,6,3,3,4,1,2)

y <- c(3,3,1,1,2,2,2,1,2,2)

polyserial(x, y, ML=TRUE)

#polyserialを用いる。

#ML=TRUEは最尤推定法、FALSEは2段階推定法

 

#########3つ以上の変数間の相関行列#########

#集計データから個票データに変換

 

individualdata <- data.frame(lapply(DanishWelfare,

                                   function(i) rep(i,DanishWelfare[,"Freq"])))

individualdata <- individualdata[,-1]

head(individualdata)

#2つのカテゴリカル変数間の相関係数を求める

polychor(individualdata$Alcohol, individualdata$Income)

 

#hectorを用いて3つ以上の変数間の相関を求める

 

#hetcorについて

#2つの変数がともにcontinuousならピアソンの積率相関係数

#ともにordinalの場合はポリコリック相関

#continuousとordinalの場合はポリシリアル相関

 

library(polycor)

ans <- hetcor(individualdata, ML=TRUE)

ans$correlations

 

#相関行列を用いて因子分析

factanal(individualdata, 1, covmat=ans$correlations)

 

#########Practice############

library(vcd)

head(Hospital)

HP <- data.frame(Hospital)

head(HP)

 

HP1 <- data.frame(lapply(HP,function(i)

  rep(i,HP[,"Freq"])))

HP1 <- HP[,-3]

head(HP1)

polychor(HP1$Visit.frequency,HP1$Length.of.stay,ML=TRUE)

plot(HP1$Visit.frequency,HP1$Length.of.stay)

 

head(PlantGrowth)

polyserial(PlantGrowth$weight,PlantGrowth$group,ML=TRUE)