Rで一般線形モデル（GLM)

##線形モデルの当てはめが出来ない

例えば不良率や生存率などのモデル化に線形モデルの仮定は成立しない。
説明変数の値によって[0,1] の区間を逸脱することからも、パラメータ推定が妥当でないこともわかる。
年齢のデータのように正の値かつ裾をひいた分布では、正規分布も当てはまらない。
そのため、非正規分布でもパラメータ推定を実行できるように拡張したのが一般化線形モデルである。

観測データの背景にあると仮定する確率分布、当てはめる関数(link関数)を用意してあげれば、R上では簡単に

##GLMの分析事例

NYのAmerican comunity surveyのデータを事例として分析します。

R

acs<- read.table("http://jaredlander.com/data/acs_ny.csv", sep=",", header=TRUE, stringsAsFactors = FALSE)

#データの変数名は以下の通り
colnames(acs)
 [1] "Acres"        "FamilyIncome" "FamilyType"   "NumBedrooms"  "NumChildren" 
 [6] "NumPeople"    "NumRooms"     "NumUnits"     "NumVehicles"  "NumWorkers"  
[11] "OwnRent"      "YearBuilt"    "HouseCosts"   "ElectricBill" "FoodStamp"   
[16] "HeatingFuel"  "Insurance"    "Language"     "income" 

#　with()で、FamiliIncomeが150000以上のデータのみを対象にする
acs$income <- with(acs, FamilyIncome >= 150000)

#使用パッケージ ggplot2 usefulをロード
library(ggplot2)
library(useful)

##ggplot()でFamiliIncomeの分布を可視化
ggplot(acs, aes( x = FamilyIncome))+
  geom_density (fill="black", alpha=0.1) +
  geom_vline (xintercept = 150000)+  #x=150000をひく
  scale_x_continuous ( label = multiple.dollar , limits = c (0, 1000000))

##lm()で線形回帰

まずは線形回帰をしてみる。

income0 <- lm (income ~ HouseCosts + NumWorkers + OwnRent+
                  NumBedrooms + FamilyType ,
                data = acs)

summary(income0) #summaryで結果表示

Call:
lm(formula = income ~ HouseCosts + NumWorkers + OwnRent + NumBedrooms + 
    FamilyType, data = acs)
Residuals:
     Min       1Q   Median       3Q      Max 
-1.06024 -0.21649 -0.10865  0.05507  1.09732 
Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)         -2.841e-01  1.069e-02 -26.580  < 2e-16 ***
HouseCosts           1.201e-04  2.161e-06  55.564  < 2e-16 ***
NumWorkers           6.026e-02  2.997e-03  20.108  < 2e-16 ***
OwnRentOutright      2.285e-01  2.895e-02   7.893 3.09e-15 ***
OwnRentRented       -3.977e-02  7.944e-03  -5.006 5.59e-07 ***
NumBedrooms          3.418e-02  2.279e-03  14.998  < 2e-16 ***
FamilyTypeMale Head  1.284e-02  1.227e-02   1.046    0.295    
FamilyTypeMarried    1.097e-01  7.076e-03  15.500  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.358 on 22737 degrees of freedom
Multiple R-squared:  0.2014,	Adjusted R-squared:  0.2011 
F-statistic:   819 on 7 and 22737 DF,  p-value: < 2.2e-16

パラメータの推定結果を表示する。

R

#係数をplotするライブラリ
library(coefplot)

coefplot(income0)

##GLM関数によるロジスティック回帰モデル

リンク関数をlogit、二項分布を仮定するロジスティック回帰モデルを推定する。

R

income1 <- glm (income ~ HouseCosts + NumWorkers + OwnRent+
                  NumBedrooms + FamilyType ,
                data = acs, family = binomial(link = "logit"))#リンク関数logit指定
summary(income1)

Call:
glm(formula = income ~ HouseCosts + NumWorkers + OwnRent + NumBedrooms + 
    FamilyType, family = binomial(link = "logit"), data = acs)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8452  -0.6246  -0.4231  -0.1743   2.9503  

Coefficients:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)         -5.738e+00  1.185e-01 -48.421   <2e-16 ***
HouseCosts           7.398e-04  1.724e-05  42.908   <2e-16 ***
NumWorkers           5.611e-01  2.588e-02  21.684   <2e-16 ***
OwnRentOutright      1.772e+00  2.075e-01   8.541   <2e-16 ***
OwnRentRented       -8.886e-01  1.002e-01  -8.872   <2e-16 ***
NumBedrooms          2.339e-01  1.683e-02  13.895   <2e-16 ***
FamilyTypeMale Head  3.336e-01  1.472e-01   2.266   0.0235 *  
FamilyTypeMarried    1.405e+00  8.704e-02  16.143   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 22808  on 22744  degrees of freedom
Residual deviance: 18073  on 22737  degrees of freedom
AIC: 18089

Number of Fisher Scoring iterations: 6

coefplot(income1)

(線形モデルとあまり変わらない・・・)

logit回帰の係数を解釈するには、逆ロジット変換を行う必要があるので

R

invlogit <- function (x) #変換関数を作成
{
  1/(1+exp(-x))
  
}
invlogit(income1$coefficients)
        (Intercept)          HouseCosts          NumWorkers     OwnRentOutright       OwnRentRented 
        0.003211572         0.500184950         0.636702036         0.854753527         0.291408659 
        NumBedrooms FamilyTypeMale Head   FamilyTypeMarried 
        0.558200010         0.582624773         0.802983719 

R

income2 <- glm (income ~ HouseCosts + NumWorkers + OwnRent+
                  NumBedrooms + FamilyType ,
                data = acs, family = binomial(link = "probit"))#リンク関数probit指定

summary(income2)

Call:
glm(formula = income ~ HouseCosts + NumWorkers + OwnRent + NumBedrooms + 
    FamilyType, family = binomial(link = "probit"), data = acs)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.9878  -0.6322  -0.4221  -0.1386   3.1324  

Coefficients:
                      Estimate Std. Error z value Pr(>|z|)    
(Intercept)         -3.235e+00  6.102e-02 -53.012   <2e-16 ***
HouseCosts           4.228e-04  9.537e-06  44.327   <2e-16 ***
NumWorkers           3.135e-01  1.429e-02  21.942   <2e-16 ***
OwnRentOutright      9.900e-01  1.190e-01   8.318   <2e-16 ***
OwnRentRented       -4.728e-01  5.048e-02  -9.365   <2e-16 ***
NumBedrooms          1.362e-01  9.692e-03  14.058   <2e-16 ***
FamilyTypeMale Head  1.661e-01  7.477e-02   2.221   0.0264 *  
FamilyTypeMarried    7.288e-01  4.342e-02  16.783   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 22808  on 22744  degrees of freedom
Residual deviance: 18038  on 22737  degrees of freedom
AIC: 18054

Number of Fisher Scoring iterations: 6

##使える関数・引数
###stringAsfacotrs - read.____ のような読み込みの時の引数

「stringsAsfactors()」は、文字列をファクタとして変換してｋるえる。
デフォルトは「TRUE」になっているが、この引数を「FALSE」にするという事は、すなわち「文字列をファクタに変換しない」という指定をする。
ファクタを使うと文字列が数値にコード化されるため、カテゴリ変数を統計モデルに導入するのに便利。
一方、ファクタは文字列のように見えるが、整数のようにふるまうため紛らわしいという点も・・・
データを読み込む際にはファクタ化せず、必要に際してファクタ化した方がエラーの原因を減らせる。