# 【3分で分かる】T法とシミュレーション実験

> 当サイト【スタビジ】の本記事では、タグチメソッドの回帰手法「T法」についてまとめていきます。重回帰分析では、説明変数から目的変数が発生していると考えています。一方でT法では、目的変数から説明変数が発生していると考えています。

- URL: https://toukei-lab.com/t%e6%b3%95
- 公開: 2022.11.29 / 更新: 2024.12.10
- サイト: スタビジ（https://toukei-lab.com）

:::talk{who="umatan"}
本記事では、タグチメソッドの回帰手法「T法」についてまとめていきます。[重回帰分析](/%e5%9b%9e%e5%b8%b0%e5%88%86%e6%9e%90)では、説明変数から目的変数が発生していると考えています。一方でT法では、目的変数から説明変数が発生していると考えています。
:::

T法とは、<mark>タグチ流の多変量解析手法の一つであり予測のための方法</mark>です。

T法について細かく見ていきましょう！

## T法とは？

用いられる、データタイプは[重回帰分析](/%e5%9b%9e%e5%b8%b0%e5%88%86%e6%9e%90)と同じです。

つまり、一つの目的変数と複数の説明変数があり、連続変量を考えています。

**重回帰分析との違い**は、<mark>データ発生の背景</mark>にあります。

重回帰分析では、説明変数から目的変数が発生していると考えています。

一方でT法では、<mark>目的変数から説明変数が発生している</mark>と考えています。\
![](/wp-content/uploads/2017/10/64830bc15a1de7807b6bcec549171a75-300x118.png)

左図：重回帰のモデル（説明変数は独立）、右図：T法のモデル

T法では、「重り」と「秤」のモデルとしてよく知られています。

一つの「重り」を複数の「秤」で測ったときにそれぞれの秤の重要性によって総合して「重り」の重さを推定する方法と考えられています。

秤の重要性は**SN比**と呼ばれます。\
![](/wp-content/uploads/2017/10/2d782027caea96e1f8a9ec4dd614af37-300x199.png)

図：重りと秤の図

## T法の解析手順

T法による解析の手順は以下の通りです。

1. 学習データの中から目的変数が中位にある$m$個のデータを単位空間とします。学習データから単位空間の平均を引くことでデータの基準化を行います。\
   \begin{eqnarray\*}\
   \bar{x}\_{j}&=&\frac{1}{m}\sum\_{i=1}^n{x^{unit}\_{ij}}, \bar{y}=\frac{1}{m}\sum\_{i=1}^n{y^{unit}\_{i}}\\\\\
   X\_{ij}&=\&x^{signal}\_{ij}-\bar{x}\_{j}, M\_{i}=y^{signal}\_{i}-\bar{y}\\\\\
   \end{eqnarray\*}

2. 比例定数$\hat{\beta_j}$とSN比$\hat{\eta_j}$を算出します。\
   \begin{eqnarray\*}\
   \hat{\beta_j}&=&\frac{\sum\_{i=1}^lX\_{ij}M_i}{r}\\\\\
   \hat{\eta_j}&=&\
   \begin{cases}\
   \frac{\frac{1}{r}(S\_{\beta_j}-V\_{e_j})}{V\_{e_j}} & (S\_{\beta_j}>V\_{e_j})\\\\\
   0 & (S\_{\beta_j} \leq V\_{e_j})\
   \end{cases} \\\\\
   r&=&{\sum\_{i=1}^lM\_{i}^2}\\\\\
   S\_{T_j}&=&{\sum\_{i=1}^lX\_{ij}^2}\\\\\
   S\_{\beta_j}&=&\frac{(\sum\_{i=1}^lX\_{ij}M_i)^2}{r}\\\\\
   S\_{e_j}&=\&S\_{T_j}-S\_{\beta_j}\\\\\
   V\_{e_j}&=&\frac{S\_{e_j}}{l-1}\\\\\
   \end{eqnarray\*}

3. 逆推定をすることで各説明変数に対する目的変数の推定値$\hat{M}_{ij}$を求めます。\
   $\displaystyle \hat{M}_{ij}=\frac{X_{ij}}{\hat{\beta_j}}$

4. 各説明変数に対する目的変数の推定値$\hat{M}_{ij}$をSN比$\hat{\eta_j}$によって重み付け総合をして基準化済みの総合推定値$\hat{M}_{i}$を算出します。\
   \begin{eqnarray\*}\
   \hat{M}\_{i}&=&\frac{1}{\sum\_{j=1}^p\hat{\eta}\_j}{\sum\_{j=1}^p{\hat{\eta}\_j\hat{M}\_{ij}}}\\\\\
   \end{eqnarray\*}

5. 単位空間の平均を加えることで基準化前の総合推定値$\hat{y}_{i}$を求めます。\
   \begin{eqnarray\*}\
   \hat{y}\_{i}&=&\frac{1}{\sum\_{j=1}^p\hat{\eta}\_j}{\sum\_{j=1}^p{\hat{\eta}\_j\hat{M}\_{ij}}}+\bar{y}\\\\\
   \end{eqnarray\*}

## T法の計算詳細

ここではT法の計算について比例項による変動とSN比の算出の仕方の説明します。

また、比例定数の二乗の推定値の求め方からもう一つのSN比が用いられることがあるのでその紹介をします。

まず、T法のモデル次は次のように式で表されます。\
\begin{eqnarray\*}\
X\_{ij}&=&\beta_jM_i+\epsilon\_{ij},\epsilon\_{ij} \sim N(0,\sigma_j^2)\\\\\
\end{eqnarray\*}

T法における比例項による変動$S_{\beta_j}$は次のように求めています。\
\begin{eqnarray\*}\
S\_{\beta_j}&=&{\sum\_{i=1}^l\hat{X}\_{ij}^2}\\\\\
&=&{\sum\_{i=1}^l{(\hat{\beta}\_jM_i)}^2}\\\\\
&=&{\sum\_{i=1}^l{\left(\frac{\sum\_{t=1}^lX\_{tj}M_t}{r}M_i\right)}^2}\\\\\
&=&{\sum\_{i=1}^l{\frac{\left({\sum\_{t=1}^lX\_{tj}M_t}\right)^2}{r^2}{M_i}^2}}\\\\\
&=&\frac{\left({\sum\_{t=1}^lX\_{tj}M_t}\right)^2}{r^2}{\sum\_{i=1}^l{{M_i}^2}}\\\\\
&=&\frac{\left({\sum\_{i=1}^lX\_{ij}M_i}\right)^2}{r}\\\\\
\end{eqnarray\*}\
続いてSN比$\eta_j$の定義は\
\begin{eqnarray\*}\
\eta_j&=&\frac{\beta_j^2}{\sigma_j^2}\
\end{eqnarray\*}\
であるので、$\hat{\eta}_j$は\
\begin{eqnarray\*}\
\hat{\eta}\_j&=&\frac{\hat{\beta}\_j^2}{\hat{\sigma}\_j^2}\
\end{eqnarray\*}\
となります。$\hat{\sigma}_j^2$は誤差分散の推定値を表すので$V_{e_j}$になります。\
そして$\hat{\beta}_j^2$は二つの計算方法があります。\
一つ目は、$\hat{\beta}_j$の二乗つまり$\beta_j$の推定値の二乗を用いる方法です。\
この場合、$\hat{\eta}_j$は\
\begin{eqnarray\*}\
\hat{\eta}\_j&=&\frac{\hat{\beta}\_j^2}{\hat{\sigma}\_j^2}\\\\\
&=&\frac{\left({\frac{\sum\_{i=1}^lX\_{ij}M_i}{r}}\right)^2}{V\_{e_j}}\\\\\
&=&\frac{{\frac{1}{r}}S\_{\beta_j}}{V\_{e_j}}\\\\\
\end{eqnarray\*}\
となります。\
二つ目は、$\beta_j^2$の推定値を用いる（推定値の二乗ではなく二乗の推定値を用いる）方法です。\
$\beta_j^2$の推定値は$\beta_j^2$と$S_{\beta_j}$の関係から次のように求めます。\
\begin{eqnarray\*}\
S\_{\beta_j}&=&\frac{\left({\sum\_{i=1}^lX\_{ij}M_i}\right)^2}{r}\\\\\
E\[S\_{\beta_j}]&=\&E\left\[\frac{\left({\sum\_{i=1}^lX\_{ij}M_i}\right)^2}{r}\right]\\\\\
&=&\frac{1}{r}E\left\[\left({\sum\_{i=1}^lX\_{ij}M_i}\right)^2\right]\\\\\
&=&\frac{1}{r}\left(E\left\[{\sum\_{i=1}^lX\_{ij}M_i}\right]^2+V\left\[{\sum\_{i=1}^lX\_{ij}M_i}\right]\right)\\\\\
\end{eqnarray\*}\
ここで、\
\begin{eqnarray\*}\
E\left\[{\sum\_{i=1}^lX\_{ij}M_i}\right]&=&{\sum\_{i=1}^l{E\[X\_{ij}]}M_i}\\\\\
&=&{\sum\_{i=1}^l(\beta_jM_i)M_i}=\beta_j{\sum\_{i=1}^lM_i^2}=\beta_jr\\\\\
V\left\[{\sum\_{i=1}^lX\_{ij}M_i}\right]&=&{\sum\_{i=1}^l{V\[X\_{ij}]}M_i^2}\\\\\
&=&{\sum\_{i=1}^l{\sigma_j^2}M_i^2}=\sigma_j^2r\\\\\
\end{eqnarray\*}\
であるので、\
\begin{eqnarray\*}\
E\[S\_{\beta_j}]&=&\frac{1}{r}\left(E\left\[{\sum\_{i=1}^lX\_{ij}M_i}\right]^2+V\left\[{\sum\_{i=1}^lX\_{ij}M_i}\right]\right)\\\\\
&=&\frac{1}{r}\left((\beta_jr)^2+(\sigma_j^2r)\right)\\\\\
&=&\beta_j^2r+\sigma_j^2\\\\\
\end{eqnarray\*}\
となります。したがって、\
\begin{eqnarray\*}\
\beta_j^2&=&\frac{1}{r}\left(E\[S\_{\beta_j}]-\sigma_j^2\right)\\\\\
\hat{\beta}\_j^2&=&\frac{1}{r}\left(S\_{\beta_j}-V\_{e_j}\right)\\\\\
\end{eqnarray\*}

この場合、$\hat{\eta}_j$は\
\begin{eqnarray\*}\
\hat{\eta}\_j&=&\frac{\hat{\beta}\_j^2}{\hat{\sigma}\_j^2}\\\\\
&=&\frac{\frac{1}{r}\left(S\_{\beta_j}-V\_{e_j}\right)}{V\_{e_j}}\\\\\
\end{eqnarray\*}

この二つの違いは式を見れば$V_{e_j}$を引くかどうかになっています。

田口玄一氏がT法を考案されたときは二つ目の方法でSN比を算出しています。

一方で、近年の論文では一つ目の方法を用いているものも多いです。

一つ目の方法は場合分けの必要がなく簡便です。

二つ目の方法は$V_{e_j}$が$S_{\beta_j}$よりも大きい場合にはその項目のSN比は0と定義するので変数選択がされます。

しかし、全てのSN比が0になる場合もありそのときには総合推定値が計算不可能になってしまいます。

結局どちらの方法も$\eta_j$の不偏推定量ではないのでどちらが良いかは分かりませんが、僕としては計算が簡便な一つ目の方法をおすすめします。

## T法を使ったシミュレーション実験

続いてT法と[重回帰分析](/%e5%9b%9e%e5%b8%b0%e5%88%86%e6%9e%90)をシミュレーション実験によって精度比較をしてみたいと思います！

まず、T法の計算について少しおさらいします。T法は個々の説明変数と目的変数とで逆単回帰を行い、変数毎の逆推定値を比例定数とSN比によって重み付け総合することで目的変数の総合推定値を算出します。総合推定値の式は次のようになりました。

\begin{eqnarray\*}\
\hat{y}\_{i}&=&\frac{1}{\sum\_{j=1}^p\hat{\eta}\_j}\sum\_{j=1}^p{\hat{\eta}\_j}\frac{X\_{ij}}{\hat{\beta_j}}+\bar{y}\\\\\
\end{eqnarray\*}

ここで、\
\begin{eqnarray\*}\
\hat{\alpha}\_{h}&=&\frac{\hat{\eta}\_h}{\hat{\beta_h}\sum\_{j=1}^p\hat{\eta}\_j}\\\\\
\end{eqnarray\*}

とおくと、\
\begin{eqnarray\*}\
\hat{y}\_{i}&=&\sum\_{j=1}^p{\hat{\alpha}\_j}{X\_{ij}}+\bar{y}\\\\\
\end{eqnarray\*}

と表記出来ます。\
よって、T法の総合推定値も重回帰モデルっぽく書くことが出来るわけです。\
そこで、まず重回帰モデルの人工データで比較してみます。\
比較手法としては重回帰、Lasso、Rdigeを取り上げます。

ただ比較において注意が必要で、厳密にはT法は学習データを単位空間と信号空間に分けて解析を行います。

信号空間は単位空間の平均によって基準化され、モデルは信号空間のデータのみで構築されます。しかし、それでは重回帰やLasso・Ridgeとはモデル構築用のデータ数が異なってしまいます。

そこで、今回はT法ではなく拡張手法であるTa法を用います。Ta法は単位空間を設定せずに全ての学習データの平均を用いて基準化する方法です。T法の拡張手法に関して詳しくは以下の記事で！

<https://toukei-lab.com/ta%e6%b3%95%e3%83%bbtb%e6%b3%95>

## 実験の設定と手順

![](/wp-content/uploads/2019/09/e37be02ae3870ffbbf95e2a91fb2d4f8_s-e1568208400539.jpg)

・比較手法：重回帰・Lasso回帰・Ridge回帰・Ta法

・評価指標：MSPE(平均二乗予測誤差)

・シミュレーション回数：100回

実験の流れは簡単にこんな感じです。

1. 人工データを生成
2. データを解析
3. 指標を算出して結果を格納
4. 1\~3を100回繰り返す
5. 結果の平均とばらつきで比較

## 人工データ生成

![](/wp-content/uploads/2018/08/PAK85_MBAdesagyou20140312_TP_V-e1557018085482-1024x575.jpg)

今回のデータの生成モデルを次に示します。

\begin{eqnarray\*}\
y_i&=&1+{\bf x}\_i^T{\bf \beta}+\epsilon_i\\\\\
{\bf \beta}&=&(\beta_1,\beta_2,0,\cdots,0)^T\in {\bf R}^p\\\\\
{\bf x_i}&\sim& N({\bf \mu},{\bf \Sigma})\\\\\
{\bf \mu}&=&{\bf 0}\_p\\\\\
{\bf \Sigma}&=&{\bf I}\_p\\\\\
\epsilon_i&\sim& N(0,1^2)\\\\\
i&=&1,2,\cdots,n\\\\\
\end{eqnarray\*}

説明変数は$p$次元であり、平均0・独立で発生させています。\
また、$\beta$は$p$次元ベクトルであり、変数1と変数2のみが目的変数へ寄与があります。\
今回は、$p$を5と100の場合で実験します。

## シミュレーション開始

![](/wp-content/uploads/2018/09/goriIMGL9703_TP_V1-e1555747456954.jpg)

シミュレーション実験をしてみましょう。\
サンプルサイズ$n=n_{train}+n_{test}$はテストデータのサンプルサイズ$n_{test}$は1000として、学習データのサンプルサイズ$n_{train}$は変化させます。

今回用いるコード（$p=5,n_{train}=20$の場合）を載せておきます。\
実行する前に入っていないパッケージはインストールしておいてください。\
ちなみにT法（Ta法）についてはMTSYSというパッケージを用いています。MTSYSには[MT法](/mt%E6%B3%95)や[RT法](/rt%e6%b3%95)も入っています。\
MTSYSを用いたTa法を実行するやり方の詳細はこちらへ！

<https://toukei-lab.com/mt%e3%82%b7%e3%82%b9%e3%83%86%e3%83%a0%e3%82%92%e5%ae%9f%e8%a3%85%e5%87%ba%e6%9d%a5%e3%82%8br%e3%83%91%e3%83%83%e3%82%b1%e3%83%bc%e3%82%b8>

```text title="Ta_method_simulation"
library(MASS)
library(mvtnorm)
library(glmnet)
library(MTSYS)
library(ggplot2)
library(tidyr)

times <- 100

n_train <- 20
n_test <- 1000
p <- 5

β0 <- 1 ##定数項
β <- c(1,-1,numeric(p-2))
β <- c(β0,β)

μ <- numeric(p)
Σ <- diag(p)

σ <- 1 ##誤差標準偏差

MSPE <- data.frame(glm=numeric(times)+NaN, lasso=numeric(times)+NaN,
                   ridge=numeric(times)+NaN, Ta=numeric(times)+NaN)

for(t in 1:times){
  X <- rmvnorm((n_train+n_test),μ,Σ)
  X <- cbind(rep(1,(n_train+n_test)),X)
  ε <- rnorm((n_train+n_test),0,σ)
  y <- X%*%β+ε
  
  data <- cbind(X,y)
  data <- as.data.frame(data)
  colnames(data) <- c("Int",paste0('x', 1:p),"y")
  
  ##学習データとテストデータを作成
  data_train <- data[1:n_train,]
  data_test <- data[-(1:n_train),]
  
  ##glm(普通の線形回帰)
  result_glm <- glm(formula=y~., data=data_train[,colnames(data_train)!="Int"], family=gaussian(link="identity"))
  
  ##glmnet(lasso)
  result_lasso <- cv.glmnet(x=as.matrix(data_train[,!(colnames(data_train) %in% c("Int","y"))]), 
                            y=data_train$y, family="gaussian",alpha=1)
  
  ##glmnet(ridge)
  result_ridge <- cv.glmnet(x=as.matrix(data_train[,!(colnames(data_train) %in% c("Int","y"))]), 
                            y=data_train$y, family="gaussian",alpha=0)
  
  ##MTSYS(Ta)
  result_Ta <- Ta(sample_data = data_train[, colnames(data_train) != "Int"],
                 subtracts_V_e = TRUE,
                 includes_transformed_data = TRUE)
  

  
  ##予測誤差を算出
  yhat_glm <- as.matrix(data_test[,colnames(data_test)!="y"])%*%as.matrix(result_glm$coefficients)
  MSPE$glm[t] <- mean((data_test$y-yhat_glm)^2)
  
  yhat_lasso <- as.matrix(data_test[,colnames(data_test)!="y"])%*%coefficients(result_lasso,s="lambda.min")
  MSPE$lasso[t] <- mean((data_test$y-yhat_lasso)^2)
  
  yhat_ridge <- as.matrix(data_test[,colnames(data_test)!="y"])%*%coefficients(result_ridge,s="lambda.min")
  MSPE$ridge[t] <- mean((data_test$y-yhat_ridge)^2)
  
  yhat_Ta <- forecasting(model = result_Ta,
                                newdata = data_test[, -which (colnames(data_test) %in% c("Int", "y"))],
                                includes_transformed_newdata = TRUE)$y_hat
  MSPE$Ta[t] <- mean((data_test$y-yhat_Ta)^2)
  
  cat("t=",t," ")
}

##ggplotで箱ひげ図を描く
g <- ggplot(gather(MSPE), aes(y=value,x=key))
g <- g + geom_boxplot()
g <- g + xlab("Method")
g <- g + ylab("MSPE")

plot(g)
```

ばらつきも含めて見たいので結果を箱ひげ図で表します。\
重回帰はglmと表記してあります。

・$p=5, n_{train}=20$の結果

![](/wp-content/uploads/2018/07/T_p5_n20.png)

やはり説明変数の数も多くないのでどの手法も同程度ですね。

・$p=5, n_{train}=1000$の結果\
学習データ数をかなり多くした結果です。

![](/wp-content/uploads/2018/07/T_p5_n1000.png)

傾向は$n_{train}=20$のときと変わりません。ここで、推定された係数を見てみましょう。

・重回帰

```text
> as.matrix(result_glm$coefficients)
[,1]
(Intercept) 0.96814597
x1 1.01957322
x2 -0.99977962
x3 0.02176648
x4 0.04343735
x5 -0.00206636
```

・Lasso

```text
> coefficients(result_lasso,s="lambda.min")
6 x 1 sparse Matrix of class "dgCMatrix"
1
(Intercept) 0.96733620
x1 0.99592767
x2 -0.97756963
x3 .
x4 0.02051056
x5 .
```

・Rdige

```text
> coefficients(result_ridge,s="lambda.min")
6 x 1 sparse Matrix of class "dgCMatrix"
1
(Intercept) 0.9675938914
x1 0.9577254944
x2 -0.9391526386
x3 0.0205584744
x4 0.0380564049
x5 -0.0004578514
```

・Ta

```text
> α_hat <- t(t(((result_Ta$eta_hat)/(result_Ta$beta_hat))/sum(result_Ta$eta_hat))) > α_hat
[,1]
x1 1.512973
x2 -1.497028
x3 0.000000
x4 0.000000
x5 0.000000
```

重回帰・Lasso・Ridgeの３手法は似たような値になっており、Lassoだけがいくつか自動的に係数が0になっています。一方で、Ta法はいくつか自動的に係数が0になっているが、変数1，2の係数が真値と比べて大きいですね。

よって、標本サイズが大きくなっても重回帰モデルのデータにおいてTa法は重回帰分析やLasso等と同じくらいの精度にはならなそうです。

・$p=100, n_{train}=10$の結果\
次に小標本高次元で実験してみます。この設定では重回帰は計算不可能なので結果から除外しています。

![](/wp-content/uploads/2018/07/T_p100_n10.png)

うーん。LassoよりTa法の方がよさそうではありますが、Ridgeと比べると微妙ですね。\
重回帰モデルの人工データでの実験なのでこんなものでしょう。

## T法を使って実データ解析

![](/wp-content/uploads/2018/09/S_8608589564326-e1566254236403-1024x575.jpg)

こんどは実データで実験してみます。\
plsパッケージにあるgasolineデータを用います。\
gasolineデータはオクタン価（１変数）とガソリンの拡散反射の近赤外線分光データ（401変数）で構成されています。\
サンプルサイズは全部で60です。\
オクタン価を近赤外線分光データで説明するモデルを構築して精度比較してみましょう。

学習データとしてランダムに30サンプルを抽出、残ったデータをテストデータとします。

1. 学習データとしてランダムに30サンプルを抽出、残ったデータをテストデータとする。
2. データを解析
3. 指標を算出して結果を格納
4. 1\~3を100回繰り返す
5. 結果の平均とばらつきで比較

実データ解析で用いたプログラムは以下の通りです。

```text title="Ta_gasoline"
library(MASS)
library(mvtnorm)
library(glmnet)
library(MTSYS)
library(ggplot2)
library(tidyr)
library(pls)

times <- 100

n_train <- 30 # < 60
n_test <- nrow(as.matrix(gasoline))-n_train
p <- ncol(as.matrix(gasoline))-1


MSPE <- data.frame(lasso=numeric(times)+NaN, ridge=numeric(times)+NaN, Ta=numeric(times)+NaN)

for(t in 1:times){
  data <- data.frame(Int = rep(1,nrow(as.matrix(gasoline))), as.matrix(gasoline)[,c((1+1):(p+1))],
                     octane = gasoline$octane)

  sample_train <- sample(1:nrow(gasoline),n_train)

  ##学習データとテストデータを作成
  data_train <- data[sample_train,]
  data_test <- data[-sample_train,]
  
  ##glmnet(lasso)
  result_lasso <- cv.glmnet(x=as.matrix(data_train[,!(colnames(data_train) %in% c("Int", "octane"))]), 
                            y=data_train$octane, family="gaussian",alpha=1)
  
  ##glmnet(ridge)
  result_ridge <- cv.glmnet(x=as.matrix(data_train[,!(colnames(data_train) %in% c("Int", "octane"))]), 
                            y=data_train$octane, family="gaussian",alpha=0)
  
  ##MTSYS(Ta)
  result_Ta <- Ta(sample_data = data_train[, colnames(data_train) != "Int"],
                  subtracts_V_e = TRUE,
                  includes_transformed_data = TRUE)
  
  ##予測誤差を算出
  yhat_lasso <- as.matrix(data_test[,colnames(data_test)!="octane"])%*%coefficients(result_lasso,s="lambda.min")
  MSPE$lasso[t] <- mean((data_test$octane-yhat_lasso)^2)
  
  yhat_ridge <- as.matrix(data_test[,colnames(data_test)!="octane"])%*%coefficients(result_ridge,s="lambda.min")
  MSPE$ridge[t] <- mean((data_test$octane-yhat_ridge)^2)
  
  yhat_Ta <- forecasting(model = result_Ta,
                         newdata = data_test[, -which (colnames(data_test) %in% c("Int", "octane"))],
                         includes_transformed_newdata = TRUE)$y_hat
  MSPE$Ta[t] <- mean((data_test$octane-yhat_Ta)^2)
  
  cat("t=",t," ")
}

##ggplotで箱ひげ図を描く
g <- ggplot(gather(MSPE), aes(y=value,x=key))
g <- g + geom_boxplot()
g <- g + xlab("Method")
g <- g + ylab("MSPE")

plot(g)
```

結果は次のようになりました。

![](/wp-content/uploads/2018/07/gasoline_Ta.png)

Ta法は全然良くないですね（笑）\
というのもこのデータは説明変数間の相関が非常に高いデータとなっています。T法（Ta法）は多重共線性が起きていても解析出来る手法ではありますが多重共線性を克服しているわけではありません。

人工データでの実験のときは説明変数間は独立であったためそこそこの精度を保てていました。

では、T法（Ta法）は使えない手法なのかというとそうとも言えません。T法は変数毎に逆単回帰を行っているだけなので逆行列の計算が必要なく非常に計算が簡便です。エクセルでも簡単に計算することが出来て、統計手法に慣れていない方でも直観的に解釈しやすい手法だと思います。

ただ、まだ適用場面が限定的であることは否めません。T法はロジスティック回帰や[ポアソン回帰](/poisson-regression)を用いるような目的変数が二値データやカウントデータへ適用する研究は不十分ですし、説明変数に質的データがあるときも適用が難しいです。

まだまだ、課題があるT法という手法ですが、解析する環境やデータの性質によっては使用すると良いかもしれません。

## T法　まとめ

本記事では、**t法の詳細やシミュレーション**について見てきました。

:::memo
**・T法（Ta法）は高次元小標本データにおいてはRidgeやLassoと同じくらいの精度が出そう。**\
**・T法のメリットは計算が簡便で解釈しやすい。**\
**・説明変数間の相関が高いと精度が悪化する恐れがある。**\
**・T法の適用場面はまだまだ限定的で研究の余地あり。**
:::

T法だけでなく、統計学や機械学習・データサイエンス、Pythonに興味のある方は以下の記事でまとめていますので是非チェックしてみてください！

<https://toukei-lab.com/statistics-freshman>

<https://toukei-lab.com/machine-learning-study>

<https://toukei-lab.com/data-scientist>

<https://toukei-lab.com/python-learn>