5.4. 非線形・交差項・二乗項#
線型モデルは制約的に感じられるかもしれませんが，二乗や自然対数ln(x)の変換を行うことで， xの非線形の効果を捉えることも可能です
5.4.1. 自然対数を用いた非線形モデルの推定#
変数の自然対数をとって推定に用いることもよく行われます。
例えば、コブダグラス型の生産関数は通常、労働\(L\)と資本\(K\)を用いてアウトプット\(Y\)が算出されるテクノロジーとして\(Y=L^aK^b\)と、定式化されます。この場合、このまま生産関数を推定しようとすると非線形モデルとなりますが、両辺とも自然対数をとって\(LnY=alnL+blnK\)と変換することで線型モデルとして扱うことができます。
また、自然対数をとると被説明変数や説明変数の単位（円、千円、時間、分など ）によらず単位の異なるデータを用いた推定結果の比較や弾性値（弾力性）の推定が容易になります。両辺が自然対数をとった\(lnY=\beta_0+\beta_1 ln X+ u\)という式ではパラメータ\(\beta_1\)は説明変数Xが1％変化した時に被説明変数yが何％変化すると解釈することができます。
傾きパラメータの解釈
被説明変数 |
説明変数 |
解釈 |
|---|---|---|
\(Y\) |
\(X\) |
\(X\)が１単位増えた時\(Y\)が\(\beta_1\)単位増える |
\(lnY\) |
\(X\) |
\(X\)が１単位増えた時\(Y\)が\(100 \times \beta_1\)%増える |
\(Y\) |
\(lnX\) |
\(X\)が1%増えた時\(Y\)が\(\beta_1 / 100\) 単位増える |
\(lnY\) |
\(lnX\) |
\(X\)が1%増えた時\(Y\)が\(\beta_1\)%増える |
5.4.2. ダミー変数（Dummy variable）#
性別や年齢層、職業カテゴリなどの連続値以外の質的データやカテゴリデータを用いて、その差が被説明変数に与える影響を調べたい際には、0または1の値を取るダミー変数を用います。
例えば、高校卒業をしているかしていないかが賃金水準に与える影響を調べたい際には、高校を卒業していれば1それ以外は0を持つダミー変数を用います。ダミー変数はこれまでに扱った連続値と同様にモデルの中に入れて推定を行うことができます。
ダミー変数を用いる際には基準カテゴリが必要となります。例えば、卒業高校に関するカテゴリ変数で、公立高校卒業、私立高校卒業、高等専門学校卒業、それ以外（高校卒業未満を含む）という4つのカテゴリを持つ変数をモデルに入れる場合、公立高校卒業者だけ1をとるダミー、私立高校卒業者だけ1をとるダミー、高等専門学校卒業者だけ1をとるダミー、それ以外（高校卒業未満を含む）だけ1をとるダミーの4つのダミーを作ることが可能ですが、実際には1つの分類を基準カテゴリ（ベース）と設定して、そのベースを除いたダミー変数のみをモデルに入れます。つまり、カテゴリの数-1このダミー変数をモデルに入れます。これは多重共線性の問題を避けるためです。
例として、賃金に経験年数が与える影響をサンプルデータを用いて見ていきます。 ソフトウェアエンジニアでの賃金の関するデータを用いて、学歴（E）と賃金(S)の関係性を分析します。 データには賃金（S）のほか、年齢（A）や経験年数（E）、学歴(E)（1=修士卒以上, 0=それ以外）、性別（G）（1=男性, 0=それ以外）が格納されています。
from statsmodels.formula.api import ols
import pandas as pd
data_path = './data/software_engineer_salary.csv'
salary_df = pd.read_csv(data_path, dtype={'A':int,'G':int,'E':int,'X':int, 'S':int})
salary_df.head(2)
| A | G | E | X | S | |
|---|---|---|---|---|---|
| 0 | 32 | 1 | 0 | 5 | 90000 |
| 1 | 28 | 0 | 0 | 3 | 125000 |
# 基本統計量を確認
salary_df['S'].describe()
count 518.000000
mean 113243.243243
std 48031.287179
min 50000.000000
25% 70000.000000
50% 90000.000000
75% 160000.000000
max 197000.000000
Name: S, dtype: float64
以下の単回帰式を求めます。
5.4.2.1. statsmodelsを用いたダミー変数の生成#
statsmodelsでは、C()でカテゴリー変数から自動的にダミー変数に割り当てます。また、Treatment()を指定することで基準カテゴリ（ベース）を指定することができます。
simple_formula = 'S ~ C(E, Treatment(0))'
# Treatmentはデフォルトでは最も若い値が基準カテゴリとなるため、
# ここでは、次の式でも同じ結果が得られます：simple_formula = 'S ~ C(EH)'
simple_res = ols(simple_formula, data = salary_df).fit()
print(simple_res.summary())
OLS Regression Results
==============================================================================
Dep. Variable: S R-squared: 0.028
Model: OLS Adj. R-squared: 0.027
Method: Least Squares F-statistic: 15.13
Date: Sat, 18 May 2024 Prob (F-statistic): 0.000114
Time: 02:19:44 Log-Likelihood: -6310.9
No. Observations: 518 AIC: 1.263e+04
Df Residuals: 516 BIC: 1.263e+04
Df Model: 1
Covariance Type: nonrobust
===========================================================================================
coef std err t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------------------
Intercept 1.156e+05 2167.488 53.327 0.000 1.11e+05 1.2e+05
C(E, Treatment(0))[T.1] -3.034e+04 7799.946 -3.889 0.000 -4.57e+04 -1.5e+04
==============================================================================
Omnibus: 5693.519 Durbin-Watson: 0.333
Prob(Omnibus): 0.000 Jarque-Bera (JB): 51.884
Skew: 0.262 Prob(JB): 5.42e-12
Kurtosis: 1.541 Cond. No. 3.77
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
上で得られた結果の中で、C(E, Treatment(0))[T.1]のcoefを確認します。
E=1（修士卒以上）のソフトウェアエンジニアは基準カテゴリである修士卒以外の従業員に比べて、賃金は$30,340平均で低いことがわかります（ただし、この分析では他の要因を考慮していないので留意が必要です）。
5.4.2.2. pandasを用いたダミー変数の生成#
次に参考までに、C()を用いないで、ダミー変数を先にpandasなどで生成してから式に入れる方法でも実行してみましょう。
ここでは、pandasのpd.get_dummies()を用いてダミー変数をテーブルデータに新たな変数として追加してから式に渡す方法を紹介します。
3つの値（1=修士卒, 2=それ以外）を持つEからダミー変数を作ります。
dummies = pd.get_dummies(salary_df['E'], prefix='E', prefix_sep='')
salary_df = salary_df.join(dummies)
# 最初の２行を表示
salary_df.head(2)
| A | G | E | X | S | E0 | E1 | |
|---|---|---|---|---|---|---|---|
| 0 | 32 | 1 | 0 | 5 | 90000 | True | False |
| 1 | 28 | 0 | 0 | 3 | 125000 | True | False |
# ダミー変数はカテゴリーの数-1の変数を式に入れ、入れなかった変数を基準カテゴリとします
simple_formula = 'S ~ E1'
simple_res = ols(simple_formula, data = salary_df).fit()
print(simple_res.summary())
OLS Regression Results
==============================================================================
Dep. Variable: S R-squared: 0.028
Model: OLS Adj. R-squared: 0.027
Method: Least Squares F-statistic: 15.13
Date: Sat, 18 May 2024 Prob (F-statistic): 0.000114
Time: 02:19:44 Log-Likelihood: -6310.9
No. Observations: 518 AIC: 1.263e+04
Df Residuals: 516 BIC: 1.263e+04
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 1.156e+05 2167.488 53.327 0.000 1.11e+05 1.2e+05
E1[T.True] -3.034e+04 7799.946 -3.889 0.000 -4.57e+04 -1.5e+04
==============================================================================
Omnibus: 5693.519 Durbin-Watson: 0.333
Prob(Omnibus): 0.000 Jarque-Bera (JB): 51.884
Skew: 0.262 Prob(JB): 5.42e-12
Kurtosis: 1.541 Cond. No. 3.77
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
statsmodelを用いた一つ前の回帰式の結果と同じになっていることが確認できます。
5.4.3. 交差項（Interaction term）#
ある変数が別の変数の値によって異なる影響を与えることを確認する際等に交差項を使います。
例えば、上の分析の例では、学歴が修士卒以上のソフトウェアエンジニアとそれ以外のソフトウェアエンジニアの間で、男性であることが賃金に与える影響の違いを捉えることができます。
上の例と同じ賃金に関するデータを用いて、交差項を見ていきます 賃金(S)と、学歴(E, 1=修士卒以上, 0=それ以外)、男性かどうか（G, 1=男性, 0=それ以外）の関係性を見ていきます。
まず、交差項を用いた分析をする前に、交差項を用いない以下の通常の重回帰分析を考えます。
formula = 'S ~ C(G, Treatment(0)) + C(E, Treatment(0))'
# formula = 'S ~ C(G) + C(E)' と書いてもここでは同じ結果が得られます。
lm = ols(formula, salary_df).fit()
print(lm.summary())
OLS Regression Results
==============================================================================
Dep. Variable: S R-squared: 0.039
Model: OLS Adj. R-squared: 0.035
Method: Least Squares F-statistic: 10.43
Date: Sat, 18 May 2024 Prob (F-statistic): 3.62e-05
Time: 02:19:44 Log-Likelihood: -6308.1
No. Observations: 518 AIC: 1.262e+04
Df Residuals: 515 BIC: 1.263e+04
Df Model: 2
Covariance Type: nonrobust
===========================================================================================
coef std err t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------------------
Intercept 1.087e+05 3614.955 30.076 0.000 1.02e+05 1.16e+05
C(G, Treatment(0))[T.1] 1.042e+04 4401.006 2.367 0.018 1769.979 1.91e+04
C(E, Treatment(0))[T.1] -2.608e+04 7971.285 -3.271 0.001 -4.17e+04 -1.04e+04
==============================================================================
Omnibus: 3922.553 Durbin-Watson: 0.339
Prob(Omnibus): 0.000 Jarque-Bera (JB): 55.259
Skew: 0.225 Prob(JB): 1.00e-12
Kurtosis: 1.464 Cond. No. 4.81
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
上の例のように交差項がない場合は、変数（男性であること）のカテゴリー間の平均値の差は、他の変数（上の例の場合学歴）に関係なく同じであると仮定しています。
交差項を用いる場合は、変数カテゴリー（上の例では男性であること）の平均値の差がが他の変数（上の例の場合は修士卒以上か）によって異なることを確認することができます。
交差項の係数の解釈では、交差項に入れた２つの変数に当てはまらい場合の平均値に比べて、追加で平均値に差をで与えると考えることができます。
次に、学歴(\(E\))と男性かどうかのダミー変数(\(G\))の交差項を考えます。
ここでは、2つのサブモデルを考えることができます。
男性でない場合(\(G=0\)): \(S_i=\beta_0+\beta_1{E_{i}}+\epsilon_i\)
男性である場合(\(G=1\)): \(S_i=(\beta_0+\beta_2)+(\beta_1+\beta_3){E_{i}}+\epsilon_i\)
パラメーターや係数の解釈例としては、以下のようになります。
\(\beta_0\)：修士卒以外（\(E=0\)）で男性ではない場合（\(G=0\)）の平均賃金
\(\beta_1\)：男性でない場合に（\(G=0\)）、修士卒（\(E=1\)）が賃金に与える効果
\((\beta_0+\beta_2)\)：修士卒以外の場合（\(E=0\)）の場合の男性であることの（\(G=1\)）が賃金に与える影響
\((\beta_1+\beta_3)\)：男性である場合（\(G=1\)）で修士卒以上である場合（\(E=1\)）が賃金に与える影響
statsmodelsでは:を用いて交差項を表現できます。C(G):C(E)で\(\beta_3E_{i1}G_{i1}\)をつくります。
interaction_formula = "S ~ C(E) + C(G) + C(E):C(G)"
interaction_ols= ols(interaction_formula, data=salary_df)
interaction_res = interaction_ols.fit()
print(interaction_res.summary())
OLS Regression Results
==============================================================================
Dep. Variable: S R-squared: 0.065
Model: OLS Adj. R-squared: 0.059
Method: Least Squares F-statistic: 11.86
Date: Sat, 18 May 2024 Prob (F-statistic): 1.61e-07
Time: 02:19:44 Log-Likelihood: -6301.0
No. Observations: 518 AIC: 1.261e+04
Df Residuals: 514 BIC: 1.263e+04
Df Model: 3
Covariance Type: nonrobust
=======================================================================================
coef std err t P>|t| [0.025 0.975]
---------------------------------------------------------------------------------------
Intercept 1.116e+05 3648.919 30.576 0.000 1.04e+05 1.19e+05
C(E)[T.1] -4.44e+04 9255.125 -4.798 0.000 -6.26e+04 -2.62e+04
C(G)[T.1] 6092.9399 4494.930 1.356 0.176 -2737.755 1.49e+04
C(E)[T.1]:C(G)[T.1] 6.624e+04 1.76e+04 3.765 0.000 3.17e+04 1.01e+05
==============================================================================
Omnibus: 4900.517 Durbin-Watson: 0.378
Prob(Omnibus): 0.000 Jarque-Bera (JB): 51.121
Skew: 0.207 Prob(JB): 7.93e-12
Kurtosis: 1.518 Cond. No. 11.0
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
得られた結果のうち交差項C(G)[T.1]:C(E)[T.1]の係数coefと、そのt値・p値を確認します。
coefは66,240で、t値が十分に大きい（p値も十分に小さい）値であることから、男性であることが賃金に与える影響は学歴（修士号以上の有無）によって効果が異なることが確認できます。
解釈としては、修士卒以外の男性以外より、修士号以上取得者は、男性であると平均で$66,240高い賃金を得ていると捉えることができます。
5.4.4. 二乗項（Quadratic term）#
賃金の上昇は若い時は経験年数が1年伸びると大きいものの、歳をとるにつれて伸びが緩やかになっていく場合もあるかもしれません。こうした変化を考慮するための１つの方法として、高次項を含む重回帰分析モデルを考えることができます。
statsmodels.from_formulaで二乗を回帰式に入れる場合はI()を直接fomula内に書くことができます。
quadr_formula = "S ~ X + I(X**2) + C(E) + C(G) + C(E):C(G)"
quadr_ols= ols(quadr_formula, data=salary_df)
quadr_res = quadr_ols.fit()
print(quadr_res.summary())
OLS Regression Results
==============================================================================
Dep. Variable: S R-squared: 0.665
Model: OLS Adj. R-squared: 0.662
Method: Least Squares F-statistic: 203.1
Date: Sat, 18 May 2024 Prob (F-statistic): 4.83e-119
Time: 02:19:44 Log-Likelihood: -6035.2
No. Observations: 518 AIC: 1.208e+04
Df Residuals: 512 BIC: 1.211e+04
Df Model: 5
Covariance Type: nonrobust
=======================================================================================
coef std err t P>|t| [0.025 0.975]
---------------------------------------------------------------------------------------
Intercept 1.528e+04 5470.497 2.792 0.005 4527.640 2.6e+04
C(E)[T.1] -1.224e+04 5745.320 -2.131 0.034 -2.35e+04 -955.646
C(G)[T.1] -3164.1318 2718.070 -1.164 0.245 -8504.075 2175.811
C(E)[T.1]:C(G)[T.1] -468.7959 1.08e+04 -0.043 0.965 -2.17e+04 2.08e+04
X 2.821e+04 1998.386 14.115 0.000 2.43e+04 3.21e+04
I(X ** 2) -995.7312 171.304 -5.813 0.000 -1332.276 -659.186
==============================================================================
Omnibus: 93.663 Durbin-Watson: 0.565
Prob(Omnibus): 0.000 Jarque-Bera (JB): 26.140
Skew: -0.264 Prob(JB): 2.11e-06
Kurtosis: 2.035 Cond. No. 337.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
得られた結果のうちI(X ** 2)が経験年数の二乗項の係数です。t値が大きくp値も小さいので、経験年数が二条項としての効果がある（経験年数が増えるほど経験年数が浅い人に比べて、経験年数1年増加に対する効果は小さくなる）と解釈できます。