[python] Logistic Regressionを用いてOdds Ratioを求める

この記事では,pythonのmoduleの一つであるstatsmodelsを用いたLogistic Regressionを行い,かつ,Odds Ratioを計算する方法を紹介する.また,Logistic Regressionまでのデータクリーニングを簡潔に行えるコードの紹介が目的である.
データクリーニングは定型作業なため,繰り返し部分を関数にして与えることによりコードの簡略化を目指すことが出来る. 

ここで用いているカスタム関数は,このブログの最後の部分に載せているコードの中に説明があるので,詳細を調べる場合はそちらを参照のこと.

動作はjupyter notebook上を想定している.
コードは,git clone https://github.com/akitoshiblog/blogsiteでlogistic_regressionの下にある.

データの準備

データの準備.最後に載せているカスタム関数に関しては,実行環境の直下に置くか,関数までのパスを追加するかを行うこと.
importするmoduleの一覧は以下の通り.

import pandas as pd
import numpy as np

from collections import Counter 
import copy
import statsmodels.api as sm
# custom functions 
import myRegression as mR

サンプルデータはstatsmodelsの中に存在するthe American National Election Studies of 1996のデータを用いる.データの詳細に関してはリンクに飛ぶと確認することが出来る.

欠損データが混入したデータセットも使用出来るように,欠損なしのデータに対しても欠損データがあるようにする.

anes96 = sm.datasets.anes96.load_pandas().data
# insert missing values
np.random.seed(130)
rnd1 = np.random.choice(anes96.index,50,replace=False) 
rnd2 = np.random.choice(anes96.index,50,replace=False) 
anes96.loc[rnd1,"age"] = np.nan
anes96.loc[rnd2,"educ"] = np.nan

データの特徴の確認

それぞれの変数が連続値かカテゴリカルなデータかを判断して入力する.
連続値ならば基本統計量を算出し,カテゴリカルなデータならばそれぞりの値に該当するデータが何件あるか出力する.

myReg = mR.myRegClass(anes96)
print(myReg.df.shape)
myReg.getValTypes(pri=True)

欠損データがどれくらい存在するかを確認する.また,そのパターンも把握出来るようにする.

myReg.checkNan()

“.getValType”のメソッドで準備した”.valTypes”のアトリビューション内に格納されている変数について情報を確認する.ここで,カテゴリー変数については,”val”キーに値が格納される.後にこの”val”キーのarrayを用いて変数をdummy変数用のカテゴリー名に変換する.

display(myReg.valTypes)

データエンジニアニングとLogist regression

the American National Election Studies of 1996のデータエンジニアニングを行う.ここでは,myRegression.pyの中で準備した機能を使うようにした.

myReg = mR.myRegClass(anes96)
myReg.getValTypes(pri=False)

# for continuous variables. 
# set the dummy categorical names. 
key = "TVnews"
bins = [0,3,6,8]
cateNames = ["seldom","moderate","frequent"]
myReg.cont2Cate(key=key,bins=bins,dummyN=cateNames,nPreKey=3)

# "eq" option specifies which place we insert equal symbol. 
key = "income"
bins = [0,8,16,24]
myReg.cont2Cate(key=key,bins=bins,dummyN=cateNames,eq="right")

# include missing values
key = "age"
bins = [19,40,60,80,100]
myReg.cont2Cate(key,bins)

key = "logpopul"
bins = [-2.5,0.09,3.09,4.70,9]
cateNames = ["high_sparse","low_sparse","low_crowded","high_crowded"]
myReg.cont2Cate(key,bins,dummyN=cateNames)

# for categorical variables. 
key = "PID"
cateNames = ["Strong Democrat","Weak Democrat","Independent-Democrat"
             ,"Independent-Indpendent", "Independent-Republican"
             ,"Weak Republican", "Strong Republican"]
myReg.getDummyVar(key = key, cateN=cateNames)

# sum up several categories into one.  
key = "selfLR"
cateNames = ["Extremely liberal","Liberal","Liberal"
             , "Moderate", "Conservative","Conservative","Extremely Conservative"]
myReg.getDummyVar(key = key, cateN=cateNames,nPreKey=4)

# includes missing values. 
key = "educ"
cateNames = ["1-8 grades","Some high school", "High school graduate"
            ,"college", "College degree", "Master's degree", "PhD"]
myReg.getDummyVar(key = key, cateN=cateNames)

最後にデータエンジニアニングが上手く出来ているか確認する.これで,欠損値がなければ問題なし.欠損値がある場合は,その変数の欠損値がダミー変数として取り込まれないことになってしまうので注意すること.
また,”cateNameSum”のアトリビューションの辞書内部には,キーに数値データをダミー変数用の名前に変更した変数名が,値には,ダミー変数のカテゴリーの名前が格納されている.

# check np.nan after replacing the name
myReg.checkNan(myReg.cateNameSum.keys())

display(myReg.df.head().T)

最後に,結果変数の選択とOdds Ratioを計算する際のベースラインを指定する辞書を用意して回すと,statsmodelsの推定結果とOdds Ratioを計算した結果が出力される.ベースラインの変数を指定しないと,説明変数を用意する際に使った一番最初の要素がベースラインとして使われる.
また,定数項はデフォルトで追加されるようになっている.

resVar = "vote"
offsetDic = {"age_c":"age_40<=_<60"}
myReg.logisticOdds(resVar,offsetDic = offsetDic)

コードの補足説明

簡単のため,ロジスティック方程式が以下のように表されるとする.


$$
z = \beta_0 + \beta_1 x_1 + \beta_2 x_2 \
p(X|\beta) = \frac{1}{1 + exp(-z)}
$$

Oddsは以下のように計算される.非常にsimpleな形になるのが特徴.


$$
\frac{ p(X|\beta)}{1 -p(X|\beta)} = exp(z)
$$

ここでx1が男性,x2が女性を表すダミー変数として,X1 = (x1,x2) = (1,0), X2 = (x1,x2) = (0,1) だとする.この時,Odds Ratioは,


$$
Odds Ratio = \frac{ p(X_1|\beta)}{1 -p(X_1|\beta)} /\frac{ p(X_2|\beta)}{1 -p(X_2|\beta)} = exp(\beta_1 - \beta_2)
$$

となる.ここでβ1もβ2も自由パラメータとすると推定が安定しないので,ベースラインに対する回帰係数は1と定めることが多い. 

ここで問題になるのは,回帰係数を1にする項をstatsmodelsのLogistic regressionの中にどう取り込むかであるが,それはoffset項として渡せば良い.offset項は線形予測子の中で係数がつかない項を指す.

以下が,"myRegression.py"の"logisticOdds"の一部で,logistic regressionをoffset項ありで行う箇所である.まずコードの上の部分でoffset項として取り扱う変数名と説明変数として取り扱う変数名を分けて保存する.
その後,offset項として取り扱う変数に関しては,一括して足し合わせて一つのpd.Seriesとして引き渡すことで実現可能.
Odds Ratioに関しては推定したパラメーターを利用して上記の式に基づいて計算することが可能である.

        # set offset variables
        cateNameSum= copy.deepcopy(self.cateNameSum)
        for k,v in cateNameSum.items(): 
            if k in offsetDic.keys():
                offsetVars.append(offsetDic[k])
                v.remove(offsetDic[k])
                expVars += v
            else:
                offsetVars.append(v.pop(0))
                expVars += v 

        offset = dfAna[offsetVars].sum(axis=1)

        # logistic regression
        mod = sm.GLM(dfAna[resVar],
                     dfAna[expVars],
                     offset=offset,
                     family=sm.families.Binomial())
        res = mod.fit()

----------雑感(`・∀・´)---------
これでLogistic Regressionのためのデータエンジニアリング,そして推定までのプロセスがシンプルに書けるようになった.
この分析を行うとわかることが,民主党の人よりも共和党の人の方が圧倒的に選挙に行きやすいということ.トランプ政権が誕生したサイレントマジョリティの存在と,民主党が選挙に勝ちにくい理由がこのデータセットからも見て取れるのは面白い! 

使用したコード

myRegression.py

#!usr/bin/python
# -*- coding: utf-8 -*-

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

import copy
from collections import Counter 

import statsmodels.api as sm



def getUnique(lis_):
    '''get the unique values and the order of given "lis_" is kept. 

    Args:
        lis_ (list) : for getting unique. 

    Return:
        list : unique list keeping the order. 
    '''
    uni_ = []
    while len(lis_):
        l = lis_.pop(0)
        uni_.append(l)
        while l in lis_:
            lis_.remove(l)
    return(uni_)
        

class myRegClass():
    def __init__(self,df):
        self.dfOri = df
        self.df = df

        self.cateNameSum = {}
        self.odds = {}


    def checkNan(self,keys=None):
        '''check Nan value patterns. 

        Args: 
            keys (list) : check the nan pattern for this variables. 
                if not given, check the all the kyes in "self.df".
            self.df : dataframe of analysis data. 

        Print:
            nan pattern.
        '''
        if keys == None:
            df_ = self.df.copy()
        else:
            keys = list(keys)
            df_ = self.df[keys].copy()
        c = Counter() 
        names = [str(s) for s in df_.columns.values]
        names = "(%s)" % ", ".join(names)
        print(names)
        
        for i in df_.index:
            c[tuple(df_.loc[i,:].isnull().replace({False:0,True:1}) ) ] +=1

        for k,v in c.items():
            print(k,v)

    def getValTypes(self, NCriteria= 20, pri = False):
        '''get types of variables, categorical or continuous. 
        types are judged by just number of kinds of data. 

        Args: 
            NCriteria (int) : more than this value is regarded as continuous. 
            pri (bool) : if True, print the types and values of each variable.
            self.df : dataframe of analysis data. 

        Sets: 
            self.valTypes (dic) : contains type, 
                and if data is categorical, possible values are included. 
                (data are sorted)
        '''
        df_ = self.df.copy() 
        dic = {}

        for k in df_.keys():
            uniq =  df_[k].unique()
            n = len(uniq)
            dic[k] = {}
            if n > NCriteria:
                dic[k]["type"] = "cont"
                if pri:
                    print("####### %s , continuous ###########" %k)
                    display(df_[k].describe().to_frame().T)
            else:
                vals = df_[k].unique()
                vals.sort()
                dic[k]["type"] = "cate"
                dic[k]["val"] = vals

                if pri:
                    print("####### %s , categorical ###########" %k)
                    cnts = df_[k].value_counts(dropna=False)
                    display(cnts[vals].to_frame().T)
        # Sets
        self.valTypes = dic

    def cont2Cate(self,key,bins,eq = "left",keyC=None,dummyN=None,nPreKey=10):
        '''continuous variables are converted into categorical variables.
        dataframe are added with dummy variables for each categorical variable.

        Args: 
            key (str) : target variable. 
            bins (list) : list of borderlines.
            eq (str) : "right" or "left" , decide "equal" is included right or left.
            keyC (str) : new categorical name. if not given, key + "_c" is applied.
            dummyN (list) : new categorical names for dummy variables.
            nPreKey (int) : number of characters from target variable for 
                prefixes of dummy variables. 
            self.df : dataframe of analysis data, which are edited after this func. 

        Sets:
            self.df : new categorical and dummy variables are added to columns. 
            self.cateNameSum (dic) : one category name and multiple names for 
                dummy variables are given.
        
        '''
        nBins = len(bins)
        if keyC == None:
            keyC = "%s_c" %key
            
        df_ = self.df

        conds = []
        cateN = []

        for i in range(nBins - 1):
            if eq == "left":
                cond = (bins[i] <= df_[key]  ) & ( df_[key] < bins[i+1])
                name = "%s_%s<=_<%s" % (key[:nPreKey],str(bins[i]),str(bins[i+1]))

            elif eq == "right":
                cond = (bins[i] < df_[key]  ) & ( df_[key] <= bins[i+1])
                name = "%s_%s<_<=%s" % (key[:nPreKey],str(bins[i]),str(bins[i+1]))

            else:
                raise Exception("'eq' takes 'left' or 'right'")


            conds.append(cond)
            cateN.append(name)

        if type(dummyN) != type(None):
            cateN = [ "%s_%s" % (key[:nPreKey], s) for s in dummyN]

        # add missing category name
        condNon = df_[key].isnull()
        if condNon.sum() > 0:
            name = "%s_miss" % (key[:nPreKey])
            conds.append(condNon)
            cateN.append(name)


        for cond, name in zip(conds, cateN):
             df_.loc[cond,keyC] = name
        dfDummy = pd.get_dummies(df_[keyC])
        df_ = pd.concat((df_,dfDummy),axis = 1 )

        # Sets
        self.df = df_
        self.cateNameSum[keyC] = getUnique(cateN)

    def getDummyVar(self,key,cateN, repV=None,keyC=None, nPreKey =10):
        '''categorical variables are converted into dummy variables.
        dataframe are added with dummy variables for each categorical variable.

        Args: 
            key (str) : target variable. 
            cateN (list) : list of categorical names for each value 
                within the target variable..
            repV (list) : list of replaced value for each "cateN". 
                same index of "cateN" and "repV" are replaced.
                if not given, values of "valTypes" are used. 
            keyC (str) : new categorical name. if not given, key + "_c" is applied.
            nPreKey (int) : number of characters from target variable for 
                prefixes of dummy variables. 
            self.df : dataframe of analysis data, which are edited after this func. 
            self.valTypes (dic) : types (and values) of each variable are contained.

        Sets:
            self.df : new categorical and dummy variables are added to columns. 
            self.cateNameSum (dic) : one category name and multiple names for 
                dummy variables are given.
        
        '''
        if keyC == None:
            keyC = "%s_c" %key
        if repV == None:
            repV = self.valTypes[key]["val"]

        df_ = self.df

        newCateN = []
        repDic = {}
        for k,v in zip(repV, cateN):
            name = "%s_%s" % ( key[:nPreKey], v)
            newCateN.append(name)
            repDic[k] = name

        # add missing category name
        condNon = df_[key].isnull()
        if condNon.sum() > 0:
            name = "%s_miss" % (key[:nPreKey])
            newCateN.append(name)
            repDic[np.nan] = name

        df_[keyC] = df_[key].replace(repDic)
        dfDummy = pd.get_dummies(df_[keyC])
        df_ = pd.concat((df_,dfDummy),axis = 1 )

        # Sets
        self.df = df_
        self.cateNameSum[keyC] = getUnique(newCateN)

    def logisticOdds(self,resVar,offsetDic = {},const=True,pri=True):
        '''logistic regression for calculating odds ratio is done. 

        Args: 
            resVar (str): result variable. Takes only 0 or 1. 
            offsetDic (dic) : determines which variable is used for baseline. 
                if not, the first category of dummy varibales is used for baseline.
            const (bool): add const term or not. 
            pri (bool) : print result or not.
            self.df : dataframe which is already added with dummy variables. 
            self.cateNameSum (dic): category and dummy variables for that category. 

        Sets:
            self.res : result of logistic regression 
            self.expVars : explanatory variables that are used in the analysis 
            self.offsetVars : explanatory variables that are used for offsets. 
            self.odds : result of odds ratios is contained here. 
        '''
        dfAna = self.df.copy()
        expVars = []
        offsetVars = []
        
        if const:
            dfAna["const"] = 1 
            expVars.append("const")

        # set offset variables
        cateNameSum= copy.deepcopy(self.cateNameSum)
        for k,v in cateNameSum.items(): 
            if k in offsetDic.keys():
                offsetVars.append(offsetDic[k])
                v.remove(offsetDic[k])
                expVars += v
            else:
                offsetVars.append(v.pop(0))
                expVars += v 

        offset = dfAna[offsetVars].sum(axis=1)

        # logistic regression
        mod = sm.GLM(dfAna[resVar],
                     dfAna[expVars],
                     offset=offset,
                     family=sm.families.Binomial())
        res = mod.fit()

        # Sets  
        self.res = res 
        self.expVars = expVars
        self.offsetVars = offsetVars

        if pri:
            print(res.summary())
            print("n#######################")
            print("##### Odds Ratios #####")
            print("#######################")

        for k in self.cateNameSum:
            self.calcOdds(k,pri=pri)
        


    def calcOdds(self,key,pri=True):
        '''calculate odds ratio using the output of logistic regression ("res")

        Args:
            key (str): target variable.
            pri (bool): print results of not.  
        '''
        self.odds[key] = {} 
        
        res = self.res
        cateNSum = copy.deepcopy(self.cateNameSum)
        cateNames = cateNSum[key]

        for v in self.offsetVars:
            if v in cateNames:
                cateNames.remove(v)
                base = v
                break
        print("n%s" % key) 
        for name in cateNames:
            OR = np.exp(res.params[name] - 1)
            cL = np.exp(res.conf_int().loc[ name,0] -1 )
            cH = np.exp(res.conf_int().loc[ name,1] - 1)
            compName = "%s / %s" % (name, base)
            self.odds[key][compName] = {}

            self.odds[key][compName]["OR"] = OR
            self.odds[key][compName]["OR_l"] = cL
            self.odds[key][compName]["OR_h"] = cH 
            if pri:
                sig = ""
                if cL > 1.0 or cH < 1.0:
                    sig = "*"
                print("%s : %2.2f%s (%2.2f, %2.2f )" %(compName, OR,sig,cL,cH) )

python
スポンサーリンク
あきとしのスクラップノート

コメント

  1. […] 以下は関係する記事.・非線形最適化関数scipy.minimizeの各手法の実装例 ・[python] Logistic Regressionを用いてOdds Ratioを求める […]

タイトルとURLをコピーしました