-R- アメリカの犯罪例に Lasso を適用

スパース推定法による統計モデリング (統計学One Point)の 1.3 にある例を実行してみた結果のメモです.

アメリカ合衆国の都市における犯罪データに関して Lasso を適用する.

まずは,Lassoの代表的なパッケージである glmnet を読み込みます.glmnet パッケージのインストール方法はこちら

> ### Lasso
> library(glmnet)
Loading required package: Matrix
Loading required package: foreach
Loaded glmnet 2.0-18


基となるアメリカ合衆国の都市における犯罪データ(crime.txt *)には,以下の6つの項目を50の都市から取得した結果が記載されています.なお,crime.txtはこちらからダウンロードできるデータです.

X1 = total overall reported crime rate per 1 million residents
         人口100万人あたりの犯罪率(目的変数 Y とする)
X2 = reported violent crime rate per 100,000 residents
         人口10万人当たりの暴力的な犯罪率 ← 今回の例では使用しない
X3 = annual police funding in $/resident
         警察の年間資金($ /住民)(説明変数 X1 とする)
X4 = % of people 25 years+ with 4 yrs. of high school
         25歳以上で,高校を卒業した人の割合(説明変数 X2 とする)
X5 = % of 16 to 19 year-olds not in highschool and not highschool graduates.
         16〜19歳のうち,高校に通っていない,もしくは卒業していない人の割合(説明変数 X3 とする)
X6 = % of 18 to 24 year-olds in college
         18〜24歳のうち,大学生の割合(説明変数 X4 とする)

X7 = % of people 25 years+ with at least 4 years of college
         25歳以上で、4年制大学を卒業した人の割合(説明変数 X5 とする)

 
このデータに対して,X1のデータを目的変数(Y)とし,X3 - X7の5つを説明変数(X1 - X5)として回帰分析を行います.


まずは,データ(crime.txt)を読み込みます.

> crime <- read.table("crime.txt")  # Load crime data


"as.matrix" を用いて行列に変換します.

> crime <- as.matrix(crime)


どのデータを使用するか,また,目的変数,説明変数の設定を行います.

> X <- crime[, 3:7]  # Explanatory variable
> y <- crime[, 1]    # Objective variable


説明変数の標準化,目的変数の中心化を行います.

> X <- scale(X)      # Standardize explanatory variables
> y <- y - mean(y)   # Centralize objective variables

ただし,glmnet パッケージでは,説明変数を標準化しなくても,初期設定で標準化され,係数の推定値は元の尺度に戻して出力されるとのことです.

Lasso推定を行います.

> # Lasso estimation
> res <- glmnet(x=X, y=y)


解パス図の描画を行います.

> # Drawing Solution path
> plot(res, xvar="lambda", label=TRUE, 
+      xlab="Logarithmic value of regularization parameter", 
+      ylab="Regression coefficient", col="black", lwd=2.5) 

f:id:HidehikoMURAO:20190723211020p:plain

 

正則化パラーメータの値を 20 に固定して,係数の推定値を求めます.

> # Fix value of regularization parameter to 20
> res1 <- glmnet(x=X, y=y, lambda=20)  
> res1$beta  # Estimated value of coefficient
5 x 1 sparse Matrix of class "dgCMatrix"
          s0
V3 133.50551
V4 -25.22804
V5  19.45576
V6   .      
V7   .      

上記の例のJupyter Notebookファイルは,GitHubStatical Modeling via Sparse Estimation.ipynbの"Apply Lasso on crime data"で見ることができます.
 
..-. --- --- - -.--- - .
* crime.txt の中身
478 184 40 74 11 31 20
494 213 32 72 11 43 18
643 347 57 70 18 16 16
341 565 31 71 11 25 19
773 327 67 72 9 29 24
603 260 25 68 8 32 15
484 325 34 68 12 24 14
546 102 33 62 13 28 11
424 38 36 69 7 25 12
548 226 31 66 9 58 15
506 137 35 60 13 21 9
819 369 30 81 4 77 36
541 109 44 66 9 37 12
491 809 32 67 11 37 16
514 29 30 65 12 35 11
371 245 16 64 10 42 14
457 118 29 64 12 21 10
437 148 36 62 7 81 27
570 387 30 59 15 31 16
432 98 23 56 15 50 15
619 608 33 46 22 24 8
357 218 35 54 14 27 13
623 254 38 54 20 22 11
547 697 44 45 26 18 8
792 827 28 57 12 23 11
799 693 35 57 9 60 18
439 448 31 61 19 14 12
867 942 39 52 17 31 10
912 1017 27 44 21 24 9
462 216 36 43 18 23 8
859 673 38 48 19 22 10
805 989 46 57 14 25 12
652 630 29 47 19 25 9
776 404 32 50 19 21 9
919 692 39 48 16 32 11
732 1517 44 49 13 31 14
657 879 33 72 13 13 22
1419 631 43 59 14 21 13
989 1375 22 49 9 46 13
821 1139 30 54 13 27 12
1740 3545 86 62 22 18 15
815 706 30 47 17 39 11
760 451 32 45 34 15 10
936 433 43 48 26 23 12
863 601 20 69 23 7 12
783 1024 55 42 23 23 11
715 457 44 49 18 30 12
1504 1441 37 57 15 35 13
1324 1022 82 72 22 15 16

940 1244 66 67 26 18 16
 
Crime data

Crime data for 50 U.S cities.
The data (X1, X2, X3, X4, X5, X6, X7) are for each city.
X1 = total overall reported crime rate per 1 million residents
X2 = reported violent crime rate per 100,000 residents
X3 = annual police funding in $/resident
X4 = % of people 25 years+ with 4 yrs. of high school
X5 = % of 16 to 19 year-olds not in highschool and not highschool graduates.
X6 = % of 18 to 24 year-olds in college

X7 = % of people 25 years+ with at least 4 years of college