第5回 回帰モデル Regression Analysis 回帰モデル Linear Regression Model 単回帰モデル、重回帰モデル simple regression, multivariate regression 非線形式の回帰分析 nonlinear regression 最小二乗法（記述統計） Least square method 回帰式の適合度のチェック（分散分析）ANOVA: analysis of variance to check fitness 最小二乗推定量の性質（推測統計） パラメータ推定量の不偏性 Unbiased estimator パラメータ推定量の確率分布とｔ検定 Distribution of the estimator and Student’s t test

回帰モデル Regression Model ある変数の値を、別の変数を用いて説明 Explain target variable by other variable(s) based on correlation among them. 目的変数 Target Y 被説明変数,従属変数 説明式を作成 Build function Explanation Variable(s) 説明変数,独立変数 変数Y,実現値yi 推計値yi=f(xi) 変数X,実現値xi Y Y Z X X

記述統計学と推測統計学 Descriptive Statistics, Inference Statistics 母集団の データ Data in Population 多数データの 数学的要約･記述 Summary of many numbers 無作為抽出 Random sampling Inference Statistics 標本集団 のデータ samples 少数データの 数学的要約･記述 Summary of smaller numbers (仮想的) 母集団 population 確率的推測･記述 inference 母数(parameter)

目的変数を[記述]する Describe the Target data using other data 目的変数 Target Y 被説明変数,従属変数 説明式を作成 Build function Explanation Variable(s) 説明変数,独立変数 変数Y,実現値yi 推計値yi=f(xi) 変数X,実現値xi 記述統計学Descriptive Sta. 関数f（X)を使いYをうまく記述 Describe Y well using function f(X) (２乗和で評価するので) 残差の二乗和を小さくしたい Minimize sum of squared error. Y 残差･誤差 X 最小二乗法 Least Square Method

最小二乗法（単回帰モデル） LSQ for simple linear model 残差平方和 SSE 未知数で微分 Differenciate SSE regarding a and b. 正規方程式　Normal Equations　上式を整理して それらを解いて

最小二乗法（重回帰モデル） LSQ for multivariate linear model 残差平方和SSE 　未知数で微分すると Differentiate respect to a,b,and c 正規方程式Normal Equations　上式を整理して 単回帰モデルと同様の連立一次方程式が出来る

最小二乗法（重回帰モデル）行列表示 LSQ in matrix form 線形モデルModel 残差 Error 残差平方和SSE 　未知数で微分すると vector differentiation 正規方程式　Normal equation 　　　　　　　　　　　　より、

行列を用いた重回帰式の計算例 lm(formula = strong ~ chemA + chemB) Coefficients: (Intercept) chemA chemB 39.728 6.433 1.750

Rによる計算例 Call: lm(formula = 強さ ~ 薬剤A) Residuals: Min 1Q Median 3Q Max -2.1111 -1.3444 -0.8111 0.8222 3.6556 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 41.478 1.128 36.761 2.86e-09 *** 薬剤A 6.433 0.874 7.361 0.000154 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.141 on 7 degrees of freedom Multiple R-squared: 0.8856, Adjusted R-squared: 0.8692 F-statistic: 54.18 on 1 and 7 DF, p-value: 0.0001545 Rによる計算例 #データフレームを作成する concrete <- data.frame( + strg=c(40.0,45.8,53.0,42.3,46.7,54.1,44.2,47.1,58.0), + chemA =c(0,1,2,0,1,2,0,1,2), + chemB =c(0,0,0,1,1,1,2,2,2) ) #データフレームから変数に代入(付置)する streng <-concrete$strg chemicalA <- concrete$chemA chemicalB <- concrete$chemB #単回帰分析の結果をkekka1に入れ，その概要を出力 kekka1 <- lm(streng~chemicalA) summary(kekka1) #重回帰分析の結果をkekka2に入れ，その概要を出力 kekka2 <- lm(streng~chemicalA+chemicalB) summary(kekka2) Call: lm(formula = 強さ ~ 薬剤A + 薬剤B) Residuals: Min 1Q Median 3Q Max -2.5611 -0.3611 0.2722 0.8222 1.9056 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 39.7278 1.0076 39.426 1.78e-08 *** 薬剤A 6.4333 0.6171 10.426 4.56e-05 *** 薬剤B 1.7500 0.6171 2.836 0.0297 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.511 on 6 degrees of freedom Multiple R-squared: 0.9511, Adjusted R-squared: 0.9348 F-statistic: 58.37 on 2 and 6 DF, p-value: 0.0001168

重共線性の問題 Multi-collinearity Problem 多数の説明変数の間に相関がある場合 High correlation among explanatory variables 目的変数への効果を一意に分離できない We cannot divide the effects of two variables 係数の推計値が安定しない（直感に反する符号を取るなど） すべての観測値が ほぼ一直線上にある この直線を含むような平面であれば、どの式を使っても当てはまりにはほとんど差はない 直線上にない場所のYの予測値には大きな差がでる Y Z X

目的変数をどの程度[記述]出来たか？ How effectively does the function describle? Xiによる説明式がない場合 なる説明式がある場合 回帰（Xiで説明 できた　y　から のずれ） yiの推計値 として、 平均値 y を 使うしかない Y Y 残差･誤差 We must use the average y, if there is no regression model. X 回帰平方和 Regressed S.S. 決定係数 Deterministic Coefficient 説明できた 平方和の割合 平均値周りの 残差（全平方和） Total S.S. 残差平方和 Error S.S.

決定係数 coefficient of determination 回帰平方和 Regressed S.S. 平均値周りの 残差（全平方和） Total S.S. 残差平方和 Error S.S. 決定係数　R squared 　説明できた平方和の割合 説明変数を増やせば1に近づく 本当はあまり意味のない変数を増やしてしまう危険性あり 自由度調整済み決定係数 Adjusted R squared Residual standard error: 1.511 on 6 degrees of freedom Multiple R-squared: 0.9511, Adjusted R-squared: 0.9348

回帰による記述（説明）力の検定 ANOVA Test of overall fitness 回帰平方和が誤差平方和に比べ、大きな意味を持っているか？Does RSS have remarkable meaning comparing to ESS? 帰無仮説：回帰平方和は誤差平方和と同程度の大きさ （回帰式は、誤差に比べて大きな説明力はない） Null Hypothesis: RSS have similar meaning with ESS. (Ratio of Two Variances originated from same population obeys to F-distribution) 対立仮説：回帰平方和は誤差平方和より大きい 　（回帰式によって誤差よりもかなり大きい部分が説明できた） Alternative Hypothesis.: RSS have larger meaning than ESS.

回帰による記述（説明）力の検定例 帰無仮説：回帰平方和は誤差平方和と同程度の大きさ （回帰式は、誤差に比べて大きな説明力はない） Rejected (棄却） Null Hypothesis: RSS have similar meaning with ESS. (Ratio of Two Variances originated from same population obeys to F-distribution) 対立仮説：回帰平方和は誤差平方和より大きい採択Approved 　（回帰式によって誤差よりもかなり大きい部分が説明できた） Alternative Hypothesis.: RSS have larger meaning than ESS. Multiple R-squared: 0.9511, Adjusted R-squared: 0.9348 F-statistic: 58.37 on 2 and 6 DF, p-value: 0.0001168

推測統計学の導入（母集団の性質の推測） Inference Statistics in Regression Analysis 各観測値は、真の値を中心にばらついていたEach observed yi distributes around yi=f(Xi). 誤差が正規分布に従うと仮定 Assume Normal distribution for error (in Population) Y Y 残差･誤差 X X 関数形　　　　を決めると、誤差の実現値　　　　　　が決まる。 それらが正規分布のもとで発生する確率を計算できる。 Once you assume f(X), you obtain the realized value of errors , then you can estimate the probability of occurrence of them.

最尤推定法としての最小二乗法 LSQ method as maximum likelihood estimation 誤差(残差)が確率分布 (正規分布)をもつと仮定 Assume normal distributed errors. 残差の実現値の同時発生確率（正規分布確率密度の積）が最大になるように、関数　　　　　を決める。 Determine f(X) such that the simultaneous probability of the observed value of errors (residuals) becomes maximum. That is given by multiplication of Normal distribution densities. Y この最大化のために最終項の残差二乗和 　　　　　　　　　　　　　　を最小化するように、 関数　　　　を定める必要がある X （記述統計学で使われてきた）最小二乗法は、 （推測統計学の考え方に従う）最尤推定法でもある。

Trial in R # data making for X x1 <- rep(1,100) x2 <- rep(seq(0,9),10) x3 <- numeric(100) for(i in 0:9){ for(j in 1:10){x3[i*10+j] <- i} } # X <- cbind(x1,x2,x3) Xt <- t(X) XtX <- Xt %*% X XtXinv <- solve(XtX) beta <- c(10,1,1) # make one sample set Y Y <- X %*% beta + rnorm(100,0,1) library(rgl) plot3d(x2,x3,Y) # regression rslt <- lm(Y ~ x2 + x3) summary(rslt) objects(rslt) rslt$coefficients rslt$residuals aicr <- AIC(rslt) devia <- aicr +rslt$rank resid <- Y - X %*% rslt$coefficients LL <- sum(log(pnorm(resid,0,1))) ess <- t(resid) %*% resid

Parameter Estimation ## estimation of beta2 # produce other values for beta2 and check residuals bet2 <-numeric(200) essq <- numeric(200) LLik <- numeric(200) for(i in 1:200) { bet2[i] <- i*0.01 resid <- Y - X %*% c(10,bet2[i],1) essq[i] <- t(resid) %*% resid LLik[i] <- sum(log(pnorm(resid,0,1))) } plot(bet2,essq) plot(bet2,LLik)

AIC 赤池の情報量基準 Akaike’s Information Criterion パラメータ数＝ｋ＋1(回帰係数と残差の分散) 最大対数尤度＝ であるから， この値が小さいほうが,無駄のない効率的なモデル式ができていることを意味している. > AIC(kekka1) [1] 42.98059 > AIC(kekka2) [1] 37.32718

パラメータ推定量の確率分布 Distribution of Parameter Estimates 各サンプルの観測値は、真の回帰式から誤差分だけずれている 誤差を含む観測値から計算したパラメータの推定値も、真の値からずれている。 Parameter estimates calculated from the sample values, each of them have error from the "true value f(Xi)” can have bias. Y Y 傾きの推定値 切片の 推定値 真の 傾きは０ 真の 切片は０ X X 各サンプルの誤差が平均0分散　　の正規分布に従うことを用い、パラメータの推定量　　が従う分布を計算すれば、真の値からかけ離れている確率をチェックできる。 Calculate distribution of parameter estimate based on population error distribution.

最小二乗推定量の確率分布 Distribution of Least-square Estimates 最小二乗推定量の行列表示 観測値の誤差構造 これらより、 Where,　　　　　　　　is not stochastic variables but fixed value. 不偏性 Unbiasedness Ε takes both negative values and positive values with same probability. Average E[ε]=0 Average of Estimates equal to true value.

最小二乗推定量の確率分布(2) Distribution of Least-square Estimates(2) If 　 obeys to Normal distribution, the linear function of ε also obeys to a Normal Distribution. パラメータの分散共分散行列 Covariance Matrix of Parameter estimates j番目のパラメータの分散： 誤差の母分散　　はわからないから、残差平方和からの分散推定値で代用 ε j番目の対角要素

最小二乗推定量の確率分布(2) Distribution of Least-square Estimates(2) Probability of Normal distributing bias Based on estimate value of variance, in stead of real variance→t distibution ｔ statistics t-distribution of d.f. of n-k Under Null Hypothesis: 　(There is no effect of j-th variable) 　　　　　　　　　　　　　　　obeys t-distribution of d.f.=n-k Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 39.7278 1.0076 39.426 1.78e-08 *** 薬剤A 6.4333 0.6171 10.426 4.56e-05 *** 薬剤B 1.7500 0.6171 2.836 0.0297 *

Produce one set of # data making for X x1 <- rep(1,100) x2 <- rep(seq(0,9),10) x3 <- numeric(100) for(i in 0:9){ for(j in 1:10){x3[i*10+j] <- i} } X <- cbind(x1,x2,x3) Xt <- t(X) XtX <- Xt %*% X XtXinv <- solve(XtX) beta <- c(10,1,1) # make one sample set Y Y <- X %*% beta + rnorm(100,0,1) # regression rslt <- lm(Y ~ x2 + x3) summary(rslt) objects(rslt) rslt$coefficients rslt$residuals aicr <- AIC(rslt) devia <- aicr +rslt$rank resid <- Y - X %*% rslt$coefficients LL <- sum(log(pnorm(resid,0,1))) ess <- t(resid) %*% resid

Monte Carlo Simulation # data making for X x1 <- rep(1,100) x2 <- rep(seq(0,9),10) x3 <- numeric(100) for(i in 0:9){ for(j in 1:10){x3[i*10+j] <- i} } # X <- cbind(x1,x2,x3) Xt <- t(X) XtX <- Xt %*% X XtXinv <- solve(XtX) beta <- c(10,1,1) # prepare the storage for the estimated parameter(s) b1 <- numeric(500) b2 <- numeric(500) b3 <- numeric(500) # Repeated Monte Carlo simulation for (i in 1:500){ Y <- X %*% beta + rnorm(100,0,1) # regression rslt <- lm(Y ~ x2 + x3) brslt <- rslt$coefficients b1[i] <- brslt[1] b2[i] <- brslt[2] b3[i] <- brslt[3] hist(b1)

回帰式のチェック plot(強さ,薬剤A) plot(kekka2$fitted.values,強さ,asp=1) あるいは，

AICに基づくステップワイズ法による変数選択 すべての変数を含むモデルをlm()で推定 その結果が入ったオブジェクトに，step()関数を適用する > step(kekka2) Start: AIC=9.79 強さ ~ 薬剤A + 薬剤B Df Sum of Sq RSS AIC <none> 13.707 9.786 - 薬剤B 1 18.375 32.082 15.440 - 薬剤A 1 248.327 262.034 34.341 Call: lm(formula = 強さ ~ 薬剤A + 薬剤B) Coefficients: (Intercept) 薬剤A 薬剤B 39.728 6.433 1.750

非線形式の回帰分析 非線形の関係式でも、置き換えによって線形化できる場合がある