Rでジニ係数を計算するときに、以下のページが便利です。
http://aoki2.si.gunma-u.ac.jp/R/Gini-index.html
しかし、これは分位数データ(各群の人数が揃っているデータ)がないときにしか使えません。
そこで、人数が揃っていないときにも累積相対度数を求めてジニ係数を計算する式を作りました。
上のページを参考に作っています。
データ例: data.csv
Area,Y,X,YX
Hokkaido,300 ,10 ,30
Hokkaido,420 ,6 ,70
Aomori,56 ,7 ,8
…
Xに人口、Yに所得、YXはY/X、つまり所得/人口をイメージしています。
Gini.index<-function( x, y, # 度数ベクトル
main="", # 図のタイトル(省略時は何も書かない)
xlab="", # x 軸の名前(省略時は何も書かない)
ylab="" # y 軸の名前(省略時は何も書かない)
)
{
stopifnot(y >= 0) # yが負でないことを確認
n <- length(y) # データの個数
ya <- cumsum(y) # yの累積度数
ys <- sum(y)
ya <- c(0, ya/ys) # yの累積相対度数
stopifnot(x >= 0) # yが負でないことを確認
xa <- cumsum(x) # xの累積度数
xs <- sum(x)
xa <- c(0, xa/xs) # xの累積相対度数
plot(xa, ya, type="l", col="blue", # これを結ぶとローレンツ曲線
main=main, xlab=xlab, ylab=ylab)
abline(v=c(0, 1), h=c(0, 1), # 外周と,
coef=c(0, 1)) # 対角線(原点を通る,傾き 1 の直線)を描く
a <- xa[-1] - xa[-length(xa)] # 台形の高さ(ベクトル)
b <- ya[-1] + ya[-length(xa)] # 上底+下底(ベクトル)
# 両ベクトルの内積はローレンツ曲線の面積の2倍
return((1-a%*%b)) # ジニ係数 1からベクトル内積を引く
}
setwd("C:/Google ドライブ/R")
d<-read.csv("data.csv", header=T)
d <- d[(d$Area=="Hokkaido"),]
d <- d[order(d$YX),] # 収入/人口(Y/X)でソート
attach(d)
Gini.index(X,Y)
さらに、地域ごとや年ごとのジニ係数一覧を求める際は、パッケージplyrのddplyを用いて以下で表として出力できます。
Gini.index <- function( x, y, # 度数ベクトル
main="", # 図のタイトル(省略時は何も書かない)
xlab="", # x 軸の名前(省略時は何も書かない)
ylab="" # y 軸の名前(省略時は何も書かない)
)
{
stopifnot(y >= 0) # yが負でないことを確認
n <- length(y) # データの個数
ya <- cumsum(y) # yの累積度数
ys <- sum(y)
ya <- c(0, ya/ys) # yの累積相対度数
stopifnot(x >= 0) # yが負でないことを確認
xa <- cumsum(x) # xの累積度数
xs <- sum(x)
xa <- c(0, xa/xs) # xの累積相対度数
# ジニ係数を求めるだけなので曲線描画は省略
a <- xa[-1] - xa[-length(xa)] # 台形の高さ(ベクトル)
b <- ya[-1] + ya[-length(xa)] # 上底+下底(ベクトル)
# 両ベクトルの内積はローレンツ曲線の面積の2倍
return((1-a%*%b)) # ジニ係数 1からベクトル内積を引く
}
setwd("C:/Google ドライブ/R")
d<-read.csv("data.csv", header=T)
d <- d[order(d$YX),] # 収入/人口(Y/X)でソート
attach(d)
library(plyr)
ddply(d, "Area", summarize, Gini=Gini.index(X, Y))