no teacher method
kagaya fumitaka
2020/1/23
fumi fumi0328 https://rpubs.com/fumi/568243 # 本サイトでは主として同志社大学の講義https://www1.doshisha.ac.jp/~mjin/R/を対象にする。
ファイル名を日本語にすると様々なエラーがでる。例えば「republish」できない。
rマークダウンについては以下参照（12/14）
ロジスティック回帰と変数選択については以下のサイト参考
Rと単回帰帰分析
- https://www1.doshisha.ac.jp/~mjin/R/Chap_14/14.html
- データセットを次のように作成する。
taikei = matrix(0,10,2)
taikei[,1]<-c(165,170,172,175,170,172,183,187,180,185)
taikei[,2]<-c(50,60,65,65,70,75,80,85,90, 95)
colnames(taikei)<- c("身長","体重")
taikei## 身長 体重
## [1,] 165 50
## [2,] 170 60
## [3,] 172 65
## [4,] 175 65
## [5,] 170 70
## [6,] 172 75
## [7,] 183 80
## [8,] 187 85
## [9,] 180 90
## [10,] 185 95plot(taikei) 
コマンド plot(taikei) で図3のような直線がない散布図を作成することができる。図3に身長を横軸、体重を縦軸にした散布図を示す。図3の直線は、散布図の点の傾向(トレンド)を表すもので、回帰直線と呼ぶ。
散布図を描画するにあたり、前記で作ったリストをデータフレーム（データベース）に変換しないといけない。
以下、データフレームへの変換
taikei.F<- data.frame(taikei) taikei.lm<- lm(体重~身長,data=taikei.F)- lm(formula, data, weights, subset, na.action)
- 各引数の意味を簡潔に次に示す。 formula は回帰式に用いる被説明変数と説明変数、定数項を用いるか用いないかなどのモデルの形式である。例えば、回帰式を y = a 0+a 1x にしたいときには formula の部分はy~x と指定し、定数項を用いない回帰式 y = a 1x にしたいときには formula の部分は y~ 1+x (あるいは y~ x -1 ) と指定する。 data は回帰分析に用いたデータセットの名前である。例えば、data = taikei、あるいは data = を省略して taikei のみでもよい。 wights は用いる説明変数に重みをつける引数である。初心者は無視してもよい。特に設定しない場合は重みを付けない。 subset はデータセットの中の一部分を用いる際に、用いる部分を明示するための添字ベクトルの引数である。指定しない場合は、全てのデータを用いる。 na.action は欠損値扱いを指定する引数である。指定がない場合は、欠損値のデータを除いたデータを用いて計算を行う。よって、na.action を指定しない場合は、予測値や残差の数はデータセットから欠損値を除いた数と等しい。 関数 lm に用いるデータ形式はデータフレームである。マトリックス形式からデータフレーム形式には次のように関数 data.frame を用いて返還する。
summary(taikei.lm)##
## Call:
## lm(formula = 体重 ~ 身長, data = taikei.F)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.158 -5.370 -2.764 6.364 9.608
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -222.1647 56.7579 -3.914 0.004454 **
## 身長 1.6809 0.3224 5.213 0.000809 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.158 on 8 degrees of freedom
## Multiple R-squared: 0.7726, Adjusted R-squared: 0.7442
## F-statistic: 27.18 on 1 and 8 DF, p-value: 0.0008091summary(takei.lm) の Residuals では残差の最小値、第1四分位数、中央値、第2四分位数、最大値が返される。次のコマンドで全ての残差を返すことができる。
residuals(taikei.lm)## 1 2 3 4 5 6 7 8
## -5.178535 -3.582877 -1.944614 -6.987219 6.417123 8.055386 -5.434165 -7.157638
## 9 10
## 9.608440 6.204098- 回帰式と回帰係数
coefficients(taikei.lm)## (Intercept) 身長
## -222.164739 1.680868- 予測値
- 回帰式で計算される値を予測値、あるいは推測値と言う。関数 summaryでは予測値は返さない。予測値を返すためにはRでは関数 predictを用いる。考察の便利のため回帰分析に用いたデータ、予測値、残差を次のように一覧表にして示す。
yosokuatai<-predict(taikei.lm)
zansa<-residuals(taikei.lm)
data.frame(taikei.F,yosokuatai,zansa)## 身長 体重 yosokuatai zansa
## 1 165 50 55.17854 -5.178535
## 2 170 60 63.58288 -3.582877
## 3 172 65 66.94461 -1.944614
## 4 175 65 71.98722 -6.987219
## 5 170 70 63.58288 6.417123
## 6 172 75 66.94461 8.055386
## 7 183 80 85.43417 -5.434165
## 8 187 85 92.15764 -7.157638
## 9 180 90 80.39156 9.608440
## 10 185 95 88.79590 6.204098plot(taikei.F)
abline(taikei.lm)
- 回帰診断図
- 回帰分析ではしばしば残差を視覚的に分析を行う。Rでは回帰分析の残差などを視覚的に考察するための環境が用意されている。回帰分析の結果を次のコマンドのように関数 plot に代入するだけで図7に示す4つの異なるグラフが返される。 関数 par は作図の環境設定関数である。用いた引数 mfrow=c(2,2) は1画面に2行2列の図を作成する指定である。回帰分析の結果を関数 plot で作成した図を回帰診断図と呼ぶ。
par(mfrow=c(1,1))#par(mfrow=c(2,2),oma=c(2,2,2,1),mar=c(5,4,3,2))
plot(taikei.lm)# plot(taikei.lm,cex=1.5,pch=21,bg="blue",col="blue")
- ① 残差とフィット値のプロット 図7（a）は残差とフィット値の散布図である。ここで言うフィット値は予測値である。図の横軸が予測値、縦軸が残差となっている。図から残差の全体像を概観することができる。この図では個体 6、8、9 の残差が相対的に大きいことがわかる。
② 残差の正規Q-Qプロット 正規Q-Qプロットはデータの正規性を考察するためにデータを視覚化する方法である。Rでは関数 qqplot を用いてQ-Qプロットを作成することができる。データが正規分布に従うと、点が直線上に並べられる。通常の回帰分析では、残差が標準正規分布に従うという仮定の下で行っている。 図7(b)は標準化した残差のQ-Qプロットである。ここで用いた例では標本データの数が少ないのでその正規性に関する議論には大きな意味がない。
③ 残差の平方根プロット 図7(c)は標準化した残差の絶対値の平和根を縦軸にし、予測値を横軸とした散布図である。この図の目的も残差の変動状況を考察することである。図から個体 6、8、9 の変動が相対的に大きいことが読みとられる。
④ Cookの距離のプロット 横軸は梃子値で、縦軸は標準化残差である。赤い点線でクックの距離0.5を示している。Cook の距離は一種の距離の測度であり、Rには関数 cooks.distance が用意されている。Cook の距離は回帰分析における影響度が大きいデータの検出などに多く用いられている。Cook の距離は全てデータ用いた場合と１つのデータを除いた後求めた回帰式による予測値を用いた場合との食い違いに関する距離の測度である。Cook の距離が大きいとそのデータが回帰式による予測値に大きく影響していることを意味する。よって、Cook の距離が大きいデータは異常値である可能性がある。Cook の距離が0.5以上であれば大きいと言われている。図7(d)では個体8の Cook の距離が比較的に大きいが0.5以下であるので個体8が異常値であるとはいえない。 診断図を作成する plot では which = n の引数を用いて、4種類の図の中の1つのみを作成することができる。ここのnは1、2、3、4でそれぞれ上記の ①、②、③、④ に対応する。例えば plot (taikei.lm,which=1) であれば残差対予測値の散布図が作成される。 回帰診断図をさらに加工したいときには、次のように個別のデータを求めて散布図を作成したほうが便利である。
- ① 残差対予測値の散布図
plot(resid(taikei.lm))#残差対予測値の散布図plot(taikei.lm$resid)
abline(h=0)
- ② Q-Qプロット
qqnorm(resid(taikei.lm))#Q-Qプロットqqnorm(taikei.lm$resid)
qqline(resid(taikei.lm))#Q-Qプロットqqline(taikei.lm$resid)
- ③ Cookの距離
plot(cooks.distance(taikei.lm))
Rと重回帰分析
weight<- c(50,60,65,65,70,75,80,85,90,95)#体重,身長,ウエスト
tall<- c(165,170,172,175,170,172,183,187,180,185)
west<- c(65,68,70,65,80,85,78,79,95,97)
taikei2<- data.frame(weight,tall,west)
taikei2## weight tall west
## 1 50 165 65
## 2 60 170 68
## 3 65 172 70
## 4 65 175 65
## 5 70 170 80
## 6 75 172 85
## 7 80 183 78
## 8 85 187 79
## 9 90 180 95
## 10 95 185 97round(cor(taikei2),4)## weight tall west
## weight 1.0000 0.8790 0.8975
## tall 0.8790 1.0000 0.5944
## west 0.8975 0.5944 1.0000pairs(taikei2,pch=21,bg="red",cex=2)
(taikei2.lm<-lm(weight~.,data=taikei2))##
## Call:
## lm(formula = weight ~ ., data = taikei2)
##
## Coefficients:
## (Intercept) tall west
## -161.6698 1.0217 0.7091- Rの初心者にとって日本語文字を使用する際に特に注意すべきなことは、日本語文字以外は全て半角文字にすることである。例えば、上記の lm 関数の中の「~」「.」「,」は全て半角である。関数の中の「.」は taikei2 データの中の、被説明変数以外の全ての変数を説明変数として用いることを示す。 返された結果 Coefficients が回帰式の係数であり、これを用いた回帰式は次のように構築する。
体重 = - 161.670 + 1.022 × 身長 + 0.709 × ウエスト 回帰式の当てはまりの良さが、どの程度であるかは回帰分析の要約情報から読み取られる。
summary(taikei2.lm)##
## Call:
## lm(formula = weight ~ ., data = taikei2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.00359 -0.56141 0.07964 1.10466 1.77918
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -161.66981 13.59310 -11.89 6.75e-06 ***
## tall 1.02172 0.08960 11.40 8.95e-06 ***
## west 0.70906 0.05729 12.38 5.17e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.6 on 7 degrees of freedom
## Multiple R-squared: 0.9901, Adjusted R-squared: 0.9872
## F-statistic: 348.7 on 2 and 7 DF, p-value: 9.782e-08- 返された決定係数 (Multiple R-Squared) は0.9901、調整済みの決定係数 (Adjusted R-squared) は0.9872である。この値は、身長のみを説明変数とした単回帰の場合の0.7726、0.7442より大きく増加している。 重回帰分析の回帰診断は、単回帰の場合と同じく関数 plot (あるいは plot.lm )を用いて、回帰診断図 (diagnostic plots) を作成し考察を行うことができる。
par(mfrow=c(2,2),oma = c(1,1,2,1),mar = c(4, 4, 2, 1))
plot(taikei2.lm,pch=21,bg=2,col=2,cex=1.5)
残差の散布図からわかるように、残差が最も大きいのは約3である。Cook の距離と残差の散布図から個体1の影響が大きいことが読み取られる。実際の問題について本格的に分析行う際には、このような個体の影響について詳細に分析行うことが必要である
① 残差対予測値の散布図
plot(resid(taikei2.lm))#残差対予測値の散布図plot(taikei.lm$resid)
abline(h=0)
- ② Q-Qプロット
qqnorm(resid(taikei2.lm))#Q-Qプロットqqnorm(taikei.lm$resid)
qqline(resid(taikei2.lm))#Q-Qプロットqqline(taikei.lm$resid)
- ③ Cookの距離
plot(cooks.distance(taikei2.lm))
Rと重回帰分析2相互作用モデル
taikei2 のデータでは、説明変数身長とウエストの間にも相関関係がある。このような説明変数間の相関関係を相互作用 (interaction) と呼ぶ。当てはまりのよいモデルを作成するためには、相互作用効果を考慮した回帰分析を行うことも可能である。説明変数が2つである場合の一般式を次に示す。 = a 0 + a 1x 1 + a 2x 2 + a 3a 1x 2 Rでは相互作用を考慮した回帰係数を簡単に求めることが可能である。次に相互作用を考慮した書式を示す。関数の中の「^2」が相互作用の指定である。
(taikei2.lm2<-lm(weight~.^2,data=taikei2))##
## Call:
## lm(formula = weight ~ .^2, data = taikei2)
##
## Coefficients:
## (Intercept) tall west tall:west
## -372.94598 2.21610 3.46816 -0.01555summary(taikei2.lm2)##
## Call:
## lm(formula = weight ~ .^2, data = taikei2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3925 -0.7506 0.1545 0.5282 1.5525
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -372.94597 85.83156 -4.345 0.00485 **
## tall 2.21610 0.48648 4.555 0.00387 **
## west 3.46816 1.11361 3.114 0.02073 *
## tall:west -0.01555 0.00627 -2.480 0.04784 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.214 on 6 degrees of freedom
## Multiple R-squared: 0.9951, Adjusted R-squared: 0.9926
## F-statistic: 405.5 on 3 and 6 DF, p-value: 2.582e-07- 決定係数は0.9951、調整済みの決定係数は 0.9926で、相互作用効果を考慮していない場合の0.9901、0.9872より増加している。調整済みの決定係数が0.9926で1に近いことからこのモデルがデータに非常によく当てはまっていると言えるであろう。次に回帰診断図を示す。
plot(taikei2.lm2,pch=21,bg=2,col=2,cex=1.5)
- ① 残差対予測値の散布図
plot(resid(taikei2.lm))#残差対予測値の散布図plot(taikei.lm$resid)
abline(h=0)
- ② Q-Qプロット
qqnorm(resid(taikei2.lm))#Q-Qプロットqqnorm(taikei.lm$resid)
qqline(resid(taikei2.lm))#Q-Qプロットqqline(taikei.lm$resid)
- ③ Cookの距離
plot(cooks.distance(taikei2.lm))
変数・モデルの選択
変数の選択は、変数を入れ換えながら回帰モデルを構築し、そのモデルを比較してより当てはまりがよいモデルを採択する。よって、変数を選択することは、モデルを選択することに等しい。 モデルの選択とは、真のモデルがあって、それに近似する複数のモデルが構築された場合、複数の近似モデルの中から真のモデルに最も近いモデルを見つけ出すことである。 調整済みの決定係数を回帰モデルの選択基準とすることもできるが、Rではモデルの選択基準として広く知られている、元統計数理研究所長赤池弘次氏が提案した AIC (Akaike’s Information Criterion) を用いて回帰モデルを選択する関数が用意されている。AIC は次の式で定義さている。
AIC = － 2 × (モデルの最大対数尤度) ＋ 2 × (モデルのパラメータ数) モデル選択の場合は、AIC の 値が小さいモデルがよいモデルであると評価する。 ここではRの中にあるデータ attitude を用いて説明することにする。データ attitude は無作為に抽出した30部門に所属している35人の事務員のアンケート回答から得られたデータである。 データセットの中の数値は各質問項目に対する好意的な反応の割合 (パーセンテージ) である。データは30行7列 (変数) のデータフレームである。その７つの変数を次に示す。
① rating： 総合評価 ② complaints： 従業員の苦情の取り扱い ③ privileges： 特別な特権は許さない ④ learning： 学習の機会 ⑤ raises： 能力基づいた昇給 ⑥ critical： 加重 ⑦ advancel： 昇進
data(attitude)
attitude[1,]## rating complaints privileges learning raises critical advance
## 1 43 51 30 39 61 92 45round(cor(attitude),2)## rating complaints privileges learning raises critical advance
## rating 1.00 0.83 0.43 0.62 0.59 0.16 0.16
## complaints 0.83 1.00 0.56 0.60 0.67 0.19 0.22
## privileges 0.43 0.56 1.00 0.49 0.45 0.15 0.34
## learning 0.62 0.60 0.49 1.00 0.64 0.12 0.53
## raises 0.59 0.67 0.45 0.64 1.00 0.38 0.57
## critical 0.16 0.19 0.15 0.12 0.38 1.00 0.28
## advance 0.16 0.22 0.34 0.53 0.57 0.28 1.00pairs(attitude,panel=panel.smooth,attitude)
- 対散布図から見られるように rating と comlaints、learning、raises との相関関係は線形的であると言えよう。 そこで、まずratingを被説明変数と残りに全ての変数を説明変数とした回帰分析の情報を求めてみる。
attitude.lm1 <- lm(rating ~ ., data = attitude)
summary(attitude.lm1 )##
## Call:
## lm(formula = rating ~ ., data = attitude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.9418 -4.3555 0.3158 5.5425 11.5990
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.78708 11.58926 0.931 0.361634
## complaints 0.61319 0.16098 3.809 0.000903 ***
## privileges -0.07305 0.13572 -0.538 0.595594
## learning 0.32033 0.16852 1.901 0.069925 .
## raises 0.08173 0.22148 0.369 0.715480
## critical 0.03838 0.14700 0.261 0.796334
## advance -0.21706 0.17821 -1.218 0.235577
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.068 on 23 degrees of freedom
## Multiple R-squared: 0.7326, Adjusted R-squared: 0.6628
## F-statistic: 10.5 on 6 and 23 DF, p-value: 1.24e-05- 返された Pr(>|t|) は complaints、learning の項目を除くと、いずれも大きい値になっている。このような場合では、変数を選択して用いる必要がる。変数の優先候補はPr(>|t|)値が小さい順である。変数の選択は、変数増加法、変数減少法などの方法がある。 変数増加法とは、Pr(>|t|) 値が最も小さい変数から順次に変数を加える方法で、最適のモデルが選択されるまで変数を付け加える操作を繰り返す。変数減少法とはPr(>|t|)値が大きい変数から順次変数を除去する方法で、最適モデルが選択されるまで、変数を削除する操作を繰り返す。Rのパッケージ stats には変数を選択する関数 step が用意されている。関数 step では、モデルの選択は AIC を基準としている。次に関数 step の使用例を示す。
attitude.lm2<-step(attitude.lm1)## Start: AIC=123.36
## rating ~ complaints + privileges + learning + raises + critical +
## advance
##
## Df Sum of Sq RSS AIC
## - critical 1 3.41 1152.4 121.45
## - raises 1 6.80 1155.8 121.54
## - privileges 1 14.47 1163.5 121.74
## - advance 1 74.11 1223.1 123.24
## <none> 1149.0 123.36
## - learning 1 180.50 1329.5 125.74
## - complaints 1 724.80 1873.8 136.04
##
## Step: AIC=121.45
## rating ~ complaints + privileges + learning + raises + advance
##
## Df Sum of Sq RSS AIC
## - raises 1 10.61 1163.0 119.73
## - privileges 1 14.16 1166.6 119.82
## - advance 1 71.27 1223.7 121.25
## <none> 1152.4 121.45
## - learning 1 177.74 1330.1 123.75
## - complaints 1 724.70 1877.1 134.09
##
## Step: AIC=119.73
## rating ~ complaints + privileges + learning + advance
##
## Df Sum of Sq RSS AIC
## - privileges 1 16.10 1179.1 118.14
## - advance 1 61.60 1224.6 119.28
## <none> 1163.0 119.73
## - learning 1 197.03 1360.0 122.42
## - complaints 1 1165.94 2328.9 138.56
##
## Step: AIC=118.14
## rating ~ complaints + learning + advance
##
## Df Sum of Sq RSS AIC
## - advance 1 75.54 1254.7 118.00
## <none> 1179.1 118.14
## - learning 1 186.12 1365.2 120.54
## - complaints 1 1259.91 2439.0 137.94
##
## Step: AIC=118
## rating ~ complaints + learning
##
## Df Sum of Sq RSS AIC
## <none> 1254.7 118.00
## - learning 1 114.73 1369.4 118.63
## - complaints 1 1370.91 2625.6 138.16- このデータの場合は Pr(>|t|)が大きいも4つの変数を除去した段階での線形モデルが、当てはまりが最も良いと判断され、変数の選択が終了した。その結果、説明変数 complaints 、learningを用いたモデルが選択された。次に回帰分析の要約を出力する。
summary(attitude.lm2)##
## Call:
## lm(formula = rating ~ complaints + learning, data = attitude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.5568 -5.7331 0.6701 6.5341 10.3610
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.8709 7.0612 1.398 0.174
## complaints 0.6435 0.1185 5.432 9.57e-06 ***
## learning 0.2112 0.1344 1.571 0.128
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.817 on 27 degrees of freedom
## Multiple R-squared: 0.708, Adjusted R-squared: 0.6864
## F-statistic: 32.74 on 2 and 27 DF, p-value: 6.058e-08- 上記の情報を用いた回帰式を次に示す。
- rating = 9.8709 + 0.6435 × complaints ＋ 0.2112 × earning
plot(attitude.lm2,pch=21,bg=2,col=2,cex=1.5)
教師なし学習を使ってデータを学習してみる-主成分分析1
-http://fistofmath.sakura.ne.jp/%E6%95%99%E5%B8%AB%E3%81%AA%E3%81%97%E5%AD%A6%E7%BF%92%E3%82%92%E4%BD%BF%E3%81%A3%E3%81%A6%E3%83%87%E3%83%BC%E3%82%BF%E3%82%92%E5%AD%A6%E7%BF%92%E3%81%97%E3%81%A6%E3%81%BF%E3%82%8B/ -https://stats.biopapyrus.jp/stats/pca/ -https://www1.doshisha.ac.jp/~mjin/R/Chap_29/29.html
library(dummies)## dummies-1.5.6 provided by Decision Patternslibrary(data.table)## Warning: package 'data.table' was built under R version 3.6.2library(tidyr)
library(ranger)## Warning: package 'ranger' was built under R version 3.6.2library(xtable)
library(nnet)## Warning: package 'nnet' was built under R version 3.6.2library(e1071)## Warning: package 'e1071' was built under R version 3.6.2library(epitools)
library(car)## Loading required package: carDatalibrary(caret)## Warning: package 'caret' was built under R version 3.6.2## Loading required package: latticelibrary(ggplot2)
library(ggthemes)## Warning: package 'ggthemes' was built under R version 3.6.2library(randomForest)## randomForest 4.6-14## Type rfNews() to see new features/changes/bug fixes.##
## Attaching package: 'randomForest'## The following object is masked from 'package:ranger':
##
## importance## The following object is masked from 'package:ggplot2':
##
## marginlibrary(rgl)
library(rattle)## Warning: package 'rattle' was built under R version 3.6.2## Rattle: A free graphical interface for data science with R.
## バージョン 5.3.0 Copyright (c) 2006-2018 Togaware Pty Ltd.
## 'rattle()' と入力して、データを多角的に分析します。##
## Attaching package: 'rattle'## The following object is masked from 'package:randomForest':
##
## importance## The following object is masked from 'package:ranger':
##
## importancelibrary(rpart.plot)## Warning: package 'rpart.plot' was built under R version 3.6.2## Loading required package: rpartlibrary(rpart)
library(epitools)
library(caret)
library(ggthemes)
library(readr)## Warning: package 'readr' was built under R version 3.6.2library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:randomForest':
##
## combine## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:data.table':
##
## between, first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(cluster)主成分分析結果のプロット-http://fistofmath.sakura.ne.jp/教師なし学習を使ってデータを学習してみる/
Data <- iris
head(Data)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosairis.data <- iris[, 1:4]
iris.name <- iris[, 5]
result_pca <- iris[, 1:4]
#主成分分析
p1 <- prcomp(iris.data, scale = TRUE)
summary(p1)## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 1.7084 0.9560 0.38309 0.14393
## Proportion of Variance 0.7296 0.2285 0.03669 0.00518
## Cumulative Proportion 0.7296 0.9581 0.99482 1.00000#主成分分析結果のプロット-http://fistofmath.sakura.ne.jp/教師なし学習を使ってデータを学習してみる/
#result_pca$xで各主成分の値を返す
plot(p1$x[, 1 : 2])
#kmeans分析
result_km_2 <- kmeans(iris.data,2)
result_km_3 <- kmeans(iris.data,3)
result_km_4 <- kmeans(iris.data, 4)
result_km_2## K-means clustering with 2 clusters of sizes 97, 53
##
## Cluster means:
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 6.301031 2.886598 4.958763 1.695876
## 2 5.005660 3.369811 1.560377 0.290566
##
## Clustering vector:
## [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1
##
## Within cluster sum of squares by cluster:
## [1] 123.79588 28.55208
## (between_SS / total_SS = 77.6 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"result_km_2$cluster## [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1#table(result_km_2,result_km_2$cluster)/50
#kmeans分析結果のプロット
col_2 <- result_km_2$cluster
col_3 <- result_km_3$cluster
col_4 <- result_km_4$cluster
plot(p1$x[, 1 : 2], col = iris[, 5])
plot(p1$x[, 1 : 2], col = col_2)
plot(p1$x[, 1 : 2], col = col_3)
plot(p1$x[, 1 : 2], col = col_4)
教師なし学習を使ってデータを学習してみる-主成分分析2
temp<-c( 50, 57, 74, 94, 112, 128, 140, 147, 150, 147, 140, 128, 112, 94, 74, 57, 57, 50, 57, 74, 94, 112, 128, 140, 147, 150, 147, 140, 128, 112, 94, 74, 74, 57, 50, 57, 74, 94, 112, 128, 140, 147, 150, 147, 140, 128, 112, 94, 94, 74, 57, 50, 57, 74, 94, 112, 128, 140, 147, 150, 147, 140, 128, 112, 112, 94, 74, 57, 50, 57, 74, 94, 112, 128, 140, 147, 150, 147, 140, 128)
okamoto<-matrix(temp, 16,5, byrow=F)
colnames(okamoto)<-c("A", "B", "C", "D", "E")
oka.pc<-prcomp(okamoto)
summary(oka.pc)## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 68.0019 40.8155 3.83644 1.30318 0.30578
## Proportion of Variance 0.7332 0.2641 0.00233 0.00027 0.00001
## Cumulative Proportion 0.7332 0.9974 0.99972 0.99999 1.00000- ソフトによっては、主成分分析の結果として固有値を返すが prcomp ではその代わりに標準偏差 (standard deviation) を返す。標準偏差は、各主成分得点の標準偏差で、固有値の正の平方根に等しい。元の変数がn個あると、非ゼロである固有値および主成分(固有ベクトル)はk () 個求まる。主成分分析の目的は、できるだけ少ない主成分に、もとの変数の情報を吸収することにある。 主成分分析を行う際には、用いる第q 番目の主成分までに、もとの変数の情報がどれだけ吸収できたかが問題となる。その判断のために、寄与率と累積寄与率に関する情報が必要となる。返された結果の中の Proportion of Variance は、各標準偏差の二乗が標準偏差の二乗の合計に占める割合(寄与率)であり、Cumulative Proportion は累積寄与率である。 返された結果から分かるように第2主成分までの累積寄与率が0.997(99,7%)であり、5次元データの情報のほとんどが第1、2主成分に縮約されている。 関数 prcomp で計算された主成分(固有ベクトル)は $rotation に、列を単位に記録する。主成分は、変数の相互関係を分析するのに用いる。丘本円周データの主成分を次に示す。
oka.pc$rotation## PC1 PC2 PC3 PC4 PC5
## A 0.3648367 6.201120e-01 0.5837693 -3.397958e-01 -0.1615786
## B 0.4805596 3.397958e-01 -0.1641399 6.201120e-01 0.4920575
## C 0.5214531 -8.326673e-17 -0.5143375 1.110223e-15 -0.6808404
## D 0.4805596 -3.397958e-01 -0.1641399 -6.201120e-01 0.4920575
## E 0.3648367 -6.201120e-01 0.5837693 3.397958e-01 -0.1615786- 主成分分析では、情報がより多く含まれている上位のいくつかの主成分を用いる。分析にはそれらの棒グラフや散布図が多く用いられている。円周データでは第2主成分までの累積寄与率が99.7%に上るので、図3に第1、2主成分の散布図のみを示す。
plot(oka.pc$rotation[,1],oka.pc$rotation[,2],type="n")
text(oka.pc$rotation[,1],oka.pc$rotation[,2],colnames(okamoto))
- さて、各個体（図2に示された16個の点）の構造を合成変数で縮約した主成分得点で再現するため、次のように第2主成分までの主成分得点の散布図を作成してみる。
plot(oka.pc$x[,1],oka.pc$x[,2],type="n")
text(oka.pc$x[,1],oka.pc$x[,2],1:16)
図4で分かるように、縦軸を軸とし180度回転すると円形が若干変形されたものの、16個の点の相対位置は図2と変わらない。 ここでは5変数 (次元) のデータを2変数に縮約している。縮約した2次元のデータは、多少歪みはあるものの、元のデータ構造をある程度再現している。
主成分分析は、変数が多いとき情報の損失を最小限に押さえながら少ない合成変数に縮約する方法であるため、元のデータ構造が100%正確に再現できないことと、再現されているのは元の変数、または個体間の相対的な関係に過ぎないことを強調しておきたい。
主成分分析を行うとき、いくつの合成変数 (主成分) を用いて分析するのが適切であるかに関しては明確な決まりがない。 分散共分散行列を用いる場合は、一般的には累積寄与率70%～80%を大まかな目安とし、累積寄与率がこれを超える主成分まで用いて分析をすることが多い。
相関行列を用いた主成分分析の場合は、固有値の値が1前後になる主成分まで用いるのが1つの目安である。 関数 prcomp には引数 scale があり、データの標準化 (データのスケールを統一) が必要なときは scale = TRUE を指定する。デフォルトには scale = FALES になっている。scale = TRUEにすると元のデータの相関行列を用いた主成分に等しい。これは関数 princomp の引数 cor = TRUE の結果と対応する。 関数 prcomp を用いた主成分分析では、関数 predict を用いて、あるデータに基づいて作成した合成変数に、新しいデータを当てることが出来る。その書き式を次に示す。
predict(object,newdata,..)
たとえば、丘本の円周データの個体1～10を用いて、合成変数を作成し、残りの11～16の個体を、求めた合成関数に当てた2次元の主成分得点の散布図を図5に示す。
oka.pc1<-prcomp(okamoto[1:10,])
oka.pc2<-predict(oka.pc1,okamoto[11:16,])
plot(oka.pc1$x[,1],oka.pc1$x[,2],type="n")
text(oka.pc1$x[,1],oka.pc1$x[,2],1:10)
plot(oka.pc2[,1],oka.pc2[,2],type="n")
text(oka.pc2[,1],oka.pc2[,2],11:16)
- 関数 prcomp の結果については、主成分と主成分得点を同一の画面上で散布図を作成する biplot 関数を用いることができる。関数 biplot を用いた結果を次に示す。
biplot(oka.pc)
- 主成分分析で求まる主成分および主成分得点の正負の符号は、固有値及び固有ベクトルを求めるときのアルゴリズムが異なると逆になる場合がある。 つまり、個体の散布図を描いたときに、異なるアルゴリズムによる結果の上下が逆になったり、左右が逆になったりすることがある。主成分分析で行う分析は、変数間、個体間の絶対的関係ではなく、相対的関係であるため、分析には問題がない。 引用文献： 丘本正(1991): 多変量解析の諸方法のモデル再現性に関する数値実験、行動計量学、Vol.18、No.2、P.47-56.
Rと因子分析(同志社大学)
因子分析 (factor analysis) は、多くの変数により記述された量的データの分析方法として、1904年にスピアーマン (Spearman) によって提案された。 因子分析で扱うデータの形式は主成分分析と基本的には同じであることから、同じ場面に利用されることが多いが、手法の開発の出発点は全く異なる。 主成分分析では、変数の間の相関関係を用いて、無相関の合成変数を求めることで多くの変数を少ない変数に縮約するが、因子分析は、変数の間の相関関係から共通因子を求めることで、多くの変数を少数個の共通因子にまとめて説明することを目的としている。 因子分析は、観測データにおける変数の間の関連成分をまとめたものを共通因子 (common factor) と呼び、他の変数と関係がなく、その変数のみ持っている成分を独自因子 (unique factor) と呼ぶ。因子分析では、観測データはお互いに関連性を持っており、共通因子と独自因子に分解できることを前提としている。 例えば、表1のような観測データが図1に示すように 個の共通因子と独自因子により構成されたとすると、式1のようなモデルで表現できる。
式の中のf ik を共通因子f kの個体iの因子得点 (factor score) と呼び、共通因子の係数を因子負荷量 (factor loading) と呼ぶ。表2に因子負荷量、表3にそれに対応する因子得点を示す。
- Rの因子分析関数 Rには最尤法による因子分析の関数 factanal が実装されている。 - factanal(x,factors,rotation,scores,…)
CATB50 = read_csv("CATB50.csv")## Warning: Missing column names filled in: 'X1' [1]## Parsed with column specification:
## cols(
## X1 = col_double(),
## A = col_double(),
## B = col_double(),
## C = col_double(),
## D = col_double(),
## E = col_double(),
## F = col_double(),
## G = col_double(),
## H = col_double(),
## I = col_double(),
## J = col_double(),
## K = col_double()
## )(CATB50.fa<-factanal(CATB50,factors=4,scores = "Bartlett"))##
## Call:
## factanal(x = CATB50, factors = 4, scores = "Bartlett")
##
## Uniquenesses:
## X1 A B C D E F G H I J K
## 0.831 0.005 0.692 0.705 0.405 0.478 0.038 0.545 0.449 0.520 0.314 0.424
##
## Loadings:
## Factor1 Factor2 Factor3 Factor4
## X1 0.396
## A 0.175 0.976
## B 0.261 0.371 0.304
## C 0.315 0.403 0.138 -0.120
## D 0.664 0.286 0.260
## E 0.594 0.300 -0.272
## F 0.199 0.960
## G 0.113 0.661
## H 0.663 0.190 0.256
## I 0.295 0.570 0.254
## J 0.803 -0.106 0.149
## K 0.302 0.199 -0.131 0.654
##
## Factor1 Factor2 Factor3 Factor4
## SS loadings 2.218 1.777 1.330 1.268
## Proportion Var 0.185 0.148 0.111 0.106
## Cumulative Var 0.185 0.333 0.444 0.549
##
## Test of the hypothesis that 4 factors are sufficient.
## The chi square statistic is 23.58 on 24 degrees of freedom.
## The p-value is 0.486head(CATB50.fa)## $converged
## [1] TRUE
##
## $loadings
##
## Loadings:
## Factor1 Factor2 Factor3 Factor4
## X1 0.396
## A 0.175 0.976
## B 0.261 0.371 0.304
## C 0.315 0.403 0.138 -0.120
## D 0.664 0.286 0.260
## E 0.594 0.300 -0.272
## F 0.199 0.960
## G 0.113 0.661
## H 0.663 0.190 0.256
## I 0.295 0.570 0.254
## J 0.803 -0.106 0.149
## K 0.302 0.199 -0.131 0.654
##
## Factor1 Factor2 Factor3 Factor4
## SS loadings 2.218 1.777 1.330 1.268
## Proportion Var 0.185 0.148 0.111 0.106
## Cumulative Var 0.185 0.333 0.444 0.549
##
## $uniquenesses
## X1 A B C D E F
## 0.83130801 0.00500000 0.69151180 0.70479962 0.40524246 0.47821252 0.03802788
## G H I J K
## 0.54452939 0.44948322 0.52048568 0.31418937 0.42426668
##
## $correlation
## X1 A B C D E
## X1 1.000000000 0.36561782 0.08117519 -0.011452402 0.10223027 0.006789583
## A 0.365617817 1.00000000 0.38731016 0.178072520 0.04565655 0.069504788
## B 0.081175189 0.38731016 1.00000000 0.147598902 0.25514632 0.083982505
## C -0.011452402 0.17807252 0.14759890 1.000000000 0.30562278 0.423540168
## D 0.102230270 0.04565655 0.25514632 0.305622776 1.00000000 0.479491265
## E 0.006789583 0.06950479 0.08398251 0.423540168 0.47949126 1.000000000
## F -0.086270342 0.13866011 0.27443137 0.445569486 0.41165626 0.397389922
## G -0.093350418 0.07191196 0.28225533 -0.009513616 0.23278005 -0.191415550
## H -0.084061765 0.05026246 0.29619540 0.301672340 0.52472778 0.371537846
## I -0.023698453 0.01012554 0.22269569 0.283381011 0.40945690 0.275550885
## J 0.024419984 -0.17335663 0.07825845 0.207462083 0.57948411 0.413192490
## K 0.086189988 -0.14729828 0.17271003 0.047035643 0.41056455 0.056194075
## F G H I J K
## X1 -0.08627034 -0.093350418 -0.08406176 -0.02369845 0.02441998 0.08618999
## A 0.13866011 0.071911956 0.05026246 0.01012554 -0.17335663 -0.14729828
## B 0.27443137 0.282255332 0.29619540 0.22269569 0.07825845 0.17271003
## C 0.44556949 -0.009513616 0.30167234 0.28338101 0.20746208 0.04703564
## D 0.41165626 0.232780046 0.52472778 0.40945690 0.57948411 0.41056455
## E 0.39738992 -0.191415550 0.37153785 0.27555089 0.41319249 0.05619408
## F 1.00000000 -0.053946618 0.32253215 0.61177399 0.24560875 0.26472690
## G -0.05394662 1.000000000 0.19190388 0.07460320 0.05454315 0.39768386
## H 0.32253215 0.191903884 1.00000000 0.23316196 0.62550439 0.40507805
## I 0.61177399 0.074603199 0.23316196 1.00000000 0.39354470 0.45571712
## J 0.24560875 0.054543153 0.62550439 0.39354470 1.00000000 0.35591526
## K 0.26472690 0.397683856 0.40507805 0.45571712 0.35591526 1.00000000
##
## $criteria
## objective counts.function counts.gradient
## 0.5681002 19.0000000 19.0000000
##
## $factors
## [1] 4(CATB50.fa)##
## Call:
## factanal(x = CATB50, factors = 4, scores = "Bartlett")
##
## Uniquenesses:
## X1 A B C D E F G H I J K
## 0.831 0.005 0.692 0.705 0.405 0.478 0.038 0.545 0.449 0.520 0.314 0.424
##
## Loadings:
## Factor1 Factor2 Factor3 Factor4
## X1 0.396
## A 0.175 0.976
## B 0.261 0.371 0.304
## C 0.315 0.403 0.138 -0.120
## D 0.664 0.286 0.260
## E 0.594 0.300 -0.272
## F 0.199 0.960
## G 0.113 0.661
## H 0.663 0.190 0.256
## I 0.295 0.570 0.254
## J 0.803 -0.106 0.149
## K 0.302 0.199 -0.131 0.654
##
## Factor1 Factor2 Factor3 Factor4
## SS loadings 2.218 1.777 1.330 1.268
## Proportion Var 0.185 0.148 0.111 0.106
## Cumulative Var 0.185 0.333 0.444 0.549
##
## Test of the hypothesis that 4 factors are sufficient.
## The chi square statistic is 23.58 on 24 degrees of freedom.
## The p-value is 0.486barplot(CATB50.fa$loadings)
- 因子負荷量の散布図(第1、2因子得点の散布図 )
plot(CATB50.fa$loadings[,1:2],type="n")
text(CATB50.fa$loadings[,1:2],colnames(CATB50))
- 因子負荷量の散布図 第3、4因子得点の散布図
plot(CATB50.fa$loadings[,3:4],type="n")
text(CATB50.fa$loadings[,3:4],colnames(CATB50))
- 因子負荷量を求めるプログラムfactor1.R - http://web.sfc.keio.ac.jp/~watanabe/adstat10.htm#1
factor1<-function(r){
n<-ncol(r)
evalue<-eigen(r)$values
evector<-eigen(r)$vectors
evalue1<-evalue[round(evalue,6)>0]
nc<-length(evalue1)
cont.evalue1<-evalue1/n
d1<-matrix(0,ncol=nc,nrow=nc)
d1[col(d1)>=row(d1)]<-1
cum.evalue<-cont.evalue1%*%d1
lamda<-diag(evalue1,nc)
evector<-evector[,1:nc]
b<-apply(evector^2,2,sum)
for ( j in 1:nc) evector[,j]<-evector[,j]/sqrt(b[j])
a<-evector%*%sqrt(lamda)
list(r=round(r,3),evalue=round(evalue1,3),
cum=round(cum.evalue,3),loading=round(a,3),
evector=round(evector,3))
}
x<-matrix(c(5,4,3,2,1,4,4,3,4,3,
5,5,4,4,3,5,5,3,4,3,
3,3,2,2,1,3,3,2,2,2,
4,3,4,2,2,2,3,2,3,1,
5,5,4,3,2,2,3,2,2,1),ncol=10,byrow=T)
out<-factor1(cor(x))
print(out)## $r
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1.000 0.839 0.688 0.563 0.468 0.343 0.563 0.612 0.559 0.280
## [2,] 0.839 1.000 0.559 0.839 0.598 0.383 0.559 0.456 0.250 0.250
## [3,] 0.688 0.559 1.000 0.563 0.802 -0.086 0.250 0.102 0.280 -0.280
## [4,] 0.563 0.839 0.563 1.000 0.869 0.514 0.688 0.408 0.280 0.280
## [5,] 0.468 0.598 0.802 0.869 1.000 0.275 0.535 0.218 0.299 0.000
## [6,] 0.343 0.383 -0.086 0.514 0.275 1.000 0.943 0.910 0.767 0.959
## [7,] 0.563 0.559 0.250 0.688 0.535 0.943 1.000 0.919 0.839 0.839
## [8,] 0.612 0.456 0.102 0.408 0.218 0.910 0.919 1.000 0.913 0.913
## [9,] 0.559 0.250 0.280 0.280 0.299 0.767 0.839 0.913 1.000 0.750
## [10,] 0.280 0.250 -0.280 0.280 0.000 0.959 0.839 0.913 0.750 1.000
##
## $evalue
## [1] 5.861 2.671 0.816 0.651
##
## $cum
## [,1] [,2] [,3] [,4]
## [1,] 0.586 0.853 0.935 1
##
## $loading
## [,1] [,2] [,3] [,4]
## [1,] -0.753 -0.368 -0.454 0.303
## [2,] -0.727 -0.478 0.071 0.487
## [3,] -0.438 -0.807 -0.304 -0.254
## [4,] -0.762 -0.469 0.443 0.059
## [5,] -0.614 -0.647 0.268 -0.364
## [6,] -0.843 0.486 0.230 -0.028
## [7,] -0.965 0.203 0.122 -0.112
## [8,] -0.892 0.413 -0.179 0.037
## [9,] -0.809 0.331 -0.368 -0.316
## [10,] -0.724 0.678 0.092 0.084
##
## $evector
## [,1] [,2] [,3] [,4]
## [1,] -0.311 -0.225 -0.503 0.376
## [2,] -0.300 -0.293 0.079 0.604
## [3,] -0.181 -0.494 -0.337 -0.315
## [4,] -0.315 -0.287 0.490 0.073
## [5,] -0.253 -0.396 0.297 -0.451
## [6,] -0.348 0.297 0.255 -0.035
## [7,] -0.399 0.124 0.135 -0.138
## [8,] -0.368 0.253 -0.199 0.046
## [9,] -0.334 0.202 -0.407 -0.392
## [10,] -0.299 0.415 0.102 0.104- 最尤法により因子を求める - http://aoki2.si.gunma-u.ac.jp/R/factanal2.html
# factanal 関数の結果を整形して表示する
factanal2 <- function( dat, # データ行列
factors=0, # 抽出する因子数
rotation=c("promax", "varimax", "none"), # 因子軸の回転法
scores=c("none", "regression", "Bartlett"), # 因子得点の算出法
verbose=TRUE) # 結果の画面表示をするかどうか
{
p <- ncol(dat) # 変数の個数
n <- nrow(dat) # 行数（欠損値を含むケースも含む）
if (is.null(colnames(dat))) colnames(dat) <- paste("Var", 1:p, sep=".") # 変数名がないときにデフォルト名をつける
if (is.null(rownames(dat))) rownames(dat) <- paste("Case", 1:n, sep="-")# 行名がないときにデフォルト名をつける
dat <- subset(dat, complete.cases(dat)) # 欠損値を持つケースを除く
rotation <- match.arg(rotation) # 引数の補完
scores <- match.arg(scores) # 引数の補完
if (factors == 0) { # 抽出因子数が指定されないときは，
factors <- max(1, floor((2*p+1-sqrt((2*p+1)^2-4*(p^2-p)))/2)) # デフォルトの因子数
}
txt <- sprintf('factanal(dat, factors=factors, rotation="%s", scores=scores)', rotation) # rotation に実際の文字列を渡すためにこのようにする
result <- eval(parse(text=txt)) # 関数呼び出し
Communality <- 1-result$uniquenesses # 共通性は，斜交回転のときには因子負荷量の二乗和ではない
result$cosmetic <- cbind(result$loadings, Communality) # 共通性を付加
if (rotation!="promax") { # 斜交回転でない場合には，
SS.loadings <- colSums(result$loadings^2) # 因子負荷量の二乗和
SS.loadings <- c(SS.loadings, sum(SS.loadings)) # 総和を加える
Proportion <- SS.loadings/p*100 # 寄与率
Cum.Prop. <- cumsum(Proportion) # 累積寄与率
Cum.Prop.[factors+1] <- NA
result$cosmetic <- rbind(result$cosmetic, SS.loadings, Proportion, Cum.Prop.)
}
if (verbose == TRUE) { # 画面表示をするとき
if (result$dof) { # モデル適合度の自由度が 0 でないとき
test <- data.frame(result$STATISTIC, result$dof, result$PVAL)
colnames(test) <- c("Chi sq.", "d.f.", "P value")
rownames(test) <- ""
cat(sprintf("H0: %i factors are sufficient.\n", factors))
print(test)
}
else { # 自由度が 0 になるとき
cat(sprintf("The degrees of freedom for the model is 0 and the fit was %g\n", result$criteria[1]))
}
cat(sprintf("\nFactor loadings(rotation:%s)\n", rotation)) # 因子負荷量
print(result$cosmetic)
if (scores != "none") {
cat(sprintf("\nFactor scores(%s)\n", scores)) # 因子得点
print(result$scores)
}
}
invisible(result) # 明示的に print しないと，何も出力しない
}
dat <- matrix(c( # 10 ケース，5 変数のデータ行列例（ファイルから読んでも良い）
-1.89, -0.02, 0.42, 1.23, -1.53,
0.06, 1.81, -0.59, -0.75, -0.12,
2.58, -0.20, -1.92, -0.49, -0.35,
0.69, -0.66, -0.77, -1.92, 0.38,
-1.05, 0.07, -0.37, -0.89, -1.62,
-2.73, 1.40, 0.18, -0.09, 0.13,
0.95, 0.84, 1.19, 1.19, 0.10,
0.93, 1.17, -1.70, 1.46, -0.25,
-0.04, -0.12, -0.34, -0.24, 1.88,
0.16, 1.03, -0.05, -0.73, 0.04
), byrow=TRUE, ncol=5)
factanal2(dat)## H0: 2 factors are sufficient.
## Chi sq. d.f. P value
## 0.2545753 1 0.6138718
##
## Factor loadings(rotation:promax)
## Factor1 Factor2 Communality
## Var.1 1.0431755 0.1481079 0.99500000
## Var.2 -0.1574652 0.2685831 0.12845163
## Var.3 -0.5203595 0.1265378 0.33586077
## Var.4 0.1579708 1.0455338 0.99500000
## Var.5 0.2153349 -0.1386064 0.08782451Rとクラスター分析(同志社大学)
我々は、物事を整理整頓する際には、機能、形状などの側面から似ているものを同じのところに集めて、片付ける。これと同じくデータについてもデータ構造の側面から似ている個体を同じのグループに仕分けることが必要である場合がある。データサイエンスにおける分類のための方法は、学習 (教師、訓練) データがある分類方法と学習データがない方法に大別される。 ここで言う学習データとは、どの個体がどのグループに属するかが既知であるデータである。グループの所属を示すデータは外的基準とも呼ばれている。学習データがある場合の分類方法は、どの個体がどのグループに属するかが既知であるデータから、分類に関するモデルを作成し、そのモデルに基づいて、グループの属性が未知であるデータを最も似ていると判断されるグループに割り当てる判別分析のような方法である。 学習データがない分類方法は、どの個体がどのグループに属するかに関する事前情報がないデータについて、グループ分けする方法で、外的基準がない分類法である。通常この種の分類方法をクラスター分析と呼ぶ。ここのクラスター (cluster) とは、花やブドウなどの房の意味で、クラスター分析とは、データの構造が似ている個体を同じの房(グループ)にまとめて、そうでないものを異なる房に集めるデータの処理方法である。 学習データがない分類方法は、どの個体がどのグループに属するかに関する事前情報がないデータについて、グループ分けする方法で、外的基準がない分類法である。通常この種の分類方法をクラスター分析と呼ぶ。ここのクラスター (cluster) とは、花やブドウなどの房の意味で、クラスター分析とは、データの構造が似ている個体を同じの房 (グループ) にまとめて、そうでないものを異なる房に集めるデータの処理方法である。 本稿では、グループの形成状態を房の形で示す樹形図を用いる階層的クラスター分析の方法と、どの個体がどのグループに属するかを示す非階層的クラスター分析方法の一種である k-means 法について紹介する。
階層的クラスター分析プロセス
seiseki<-matrix(c(89,90,67,46,50, 57,70,80,85,90,80,90,35,40,50, 40,60,50,45,55,78,85,45,55,60, 55,65,80,75,85,90,85,88,92,95),7,5,byrow = TRUE)
colnames(seiseki)<- c("算数","理科","国語","英語","社会")
rownames(seiseki)<- c("田中","佐藤","鈴木","本田","川端","吉野","斉藤")
seiseki.d<- dist(seiseki)#このデータのユークリッド距離を次に示す。
round(seiseki.d)## 田中 佐藤 鈴木 本田 川端 吉野
## 佐藤 69
## 鈴木 34 81
## 本田 60 64 53
## 川端 28 61 21 47
## 吉野 63 12 76 54 56
## 斉藤 68 38 88 92 68 46階層的クラスター分析樹形図の切断とコーフエン相関係数(同志社大学)
- https://www1.doshisha.ac.jp/~mjin/R/Chap_29/29.html
- 関数 hclust パッケージ stats には階層的クラスター分析の関数 hclust がる。その書き式を次に示す。 hclust(d, method = “complete”, …)
sei.d<-dist(seiseki)
(sei.hc<-hclust(sei.d))##
## Call:
## hclust(d = sei.d)
##
## Cluster method : complete
## Distance : euclidean
## Number of objects: 7summary(sei.hc)## Length Class Mode
## merge 12 -none- numeric
## height 6 -none- numeric
## order 7 -none- numeric
## labels 7 -none- character
## method 1 -none- character
## call 2 -none- call
## dist.method 1 -none- charactersei.hc$merge#merge は、クラスター形成の履歴のマトリックスである。## [,1] [,2]
## [1,] -2 -6
## [2,] -3 -5
## [3,] -1 2
## [4,] -7 1
## [5,] -4 3
## [6,] 4 5sei.hc$height#この値は、merge の結果と対応する。例えば、個体2 (佐藤) と個体6 (吉野) の距離 (枝の長さ) ## [1] 12.40967 21.30728 33.77869 45.58509 60.13319 91.53142sei.hc$order## [1] 7 2 6 4 1 3 5- 樹形図
par(mfrow=c(2,2))
plot(sei.hc,main="最遠隣法")
plot(sei.hc,hang=-1,main="最遠隣法")
s.hc2<-hclust(sei.d,method="centroid")
plot(s.hc2,hang=-1,main="重心法")
s.hc3<-hclust(sei.d,method="ward.D")
plot(s.hc3,hang=-1,main="ウォード法")
- コーフエン相関係数-あやめのクラスター分析
- https://www1.doshisha.ac.jp/~mjin/R/Chap_28/28.html
iris2<-iris[51:150,1:4]
iris2.hc<-hclust(dist(iris2),method="ward.D")
(iris2.cl<-cutree(iris2.hc,k=2))## 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110
## 1 1 1 1 1 1 1 1 1 1 2 1 2 2 2 2 1 2 2 2
## 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
## 2 2 2 1 1 2 2 2 2 1 2 1 2 1 2 2 1 1 2 2
## 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
## 2 2 2 1 2 2 2 2 1 2 2 2 1 2 2 2 1 2 2 1- クラスター分類の結果の精度を確認するため、既知のクラス情報と分類結果のクロス表を次に示す。
iris2.lab<-c(rep(1,50),rep(2,50))
table(iris2.lab, iris2.cl)## iris2.cl
## iris2.lab 1 2
## 1 50 0
## 2 14 36- コーフエン行列を返す関数 cophenetic
iris2.cop<-cophenetic(iris2.hc)- コーフエン相関係数は、距離の行列と用いた方法で生成したコーフエン行列とのピアソン相関係数である。次にユークリッド距離とウォード法によるコーフエン相関係数を求める例を示す。
cor(iris2.cop,dist(iris2))## [1] 0.6184636- コーフエン相関係数と正解率 このコーフエン相関係数の値が大きいほど、距離行列と用いた方法のコーフエン行列との歪みが小さいと、ある側面から言える。しかし、歪みが小さいことと分類の結果がより妥当であることは等価ではない。 ここでは、データ iris2 のユークリッド (Euclidien) 距離とキャンベラ (canberra) 距離を用いたコーフエン相関係数と正しく分類された比率(正解率)を表1に示す。 正解率は、関数 cutree を用いて各個体のクラスの属性を求め、既知のクラス属性情報とのクロス表から次のように求めた。
iris2.ta<-table(iris2.lab,iris2.cl)
sum(diag(iris2.ta))/100## [1] 0.86関数 kmeans と例
iris2.km$cluster にクラスターの分類結果が記録されている。次のコマンドで、クラスターの精度を確認することができる。
iris2.km<-kmeans(iris2,2)
table(iris2.lab,iris2.km$cluster)/50##
## iris2.lab 1 2
## 1 0.96 0.04
## 2 0.28 0.72モデルに基づいたクラスター分析（同志社大学）
-https://www1.doshisha.ac.jp/~mjin/R/Chap_29/29.html - モデルに基づいたクラスター分析は、通常モデルに基づいたクラスタリング (Model-Based Clustering)、混合分布によるクラスター分析、潜在クラスター分析とも呼ばれている。 モデルに基づいたクラスター分析は、観測データが異なる分布の混合分布であると仮定し、個体が属するクラスのラベルをも隠れ変数として推定する。混合分布のパラメータおよびクラスのラベルの推定は EM (expectation maximization) アルゴリズムが多く用いられている[5]。 パッケージ mclust には関数 mclustBIC、hc、hclust などモデルに基づいたクラスタリングに関連する豊富な関数が用意されている。 関数 mclustBIC は、最大尤度推測法を用いたEMアルゴリズムでパラメータを推定する際に必要となる、ガウス混合分布モデルの情報量規準 BIC (Bayesian Information Criterion) 値を求める。通常情報量規準は、値が小さいモデルが真のモデルに、より近似していると判断するが、パッケージ mclust では、BIC の値が大きいモデルが真のモデルに、より近似すると判断する。つまり、情報量規準の式の正負が通常の式とは逆である。 関数 hc は、モデルに基づいた階層的クラスタリングを行う。関数 hc では、混合分布の型(球、楕円球)、体積、形と軸の向きが同じであるかどうかを引数 modelName =" " で指定するようになっている。多変量の場合は次の4種類から１つを選択できる。
library(mclust)## Warning: package 'mclust' was built under R version 3.6.2## Package 'mclust' version 5.4.5
## Type 'citation("mclust")' for citing this R package in publications.plot(mclustBIC(iris2))
- 図2では、折れ線Gのクラスター数2,3のBIC値が最も大きい。GはモデルEEE(楕円球型、同体積・形・向き)である。モデルEEEに基づいた階層的クラスタリングは次のように求める。
mhc<- hc(modelName = "EEE", data = iris2)- さらに、混合分布EEEモデルに基づいて求めた結果を、次のように関数 hclass にクラスターの数を指定することで、個体が属するクラスを決めることができる。
cl<- hclass(mhc,2)- クラスタリングの精度を確認するため、既知のクラスのラベルと推測結果とのクロス表を次に示す。
table(iris2.lab,cl)## cl
## iris2.lab 1 2
## 1 49 1
## 2 5 45- 関数 clPairs を次のように用いることにより、推定されたクラスターの対散布図を作成することができる。
clPairs(iris2,cl=cl)
Rによるニューラルネットワーク
- https://www1.doshisha.ac.jp/~mjin/R/Chap_23/23.html
- パッケージ nnet における関数 predict の書き式を次に示す。 - predict(学習したモデル, データ,type="")
- 引数 type を “raw” に指定すると当てはめの値を返し、“class” -に指定すると対応するグループ(判別分析のグループ)に関する値を返す。
- 実際のデータを用いて関数 nnet の使用方法について説明する。ここでも、iris
- データを用い、まず判別分析の場合と同じく次のように学習用のデータ (train.data) とテスト用のデータ (test.data) を作成する。
even.n<-2*(1:75) train.data<-iris[even.n,] test.data<-iris[-even.n,] Iris.lab<-factor(c(rep(“S”,25),rep(“C”,25),rep(“V”,25))) train.data[,5]<-Iris.lab test.data[,5]<-Iris.lab
even.n<-2*(1:75)
train.data<-iris[even.n,]
test.data<-iris[-even.n,]
Iris.lab<-factor(c(rep("S",25),rep("C",25),rep("V",25)))
train.data[,5]<-Iris.lab
test.data[,5]<-Iris.lab- 作成した学習データを用いて次のように関数 nnet によるニューラルネットワークモデルを作成し、関数 predict を用いてテストデータについて予測・判別を行なう。
library(nnet)
iris.nnet<- nnet(Species~ .,size=3,decay=0.1,data=train.data)## # weights: 27
## initial value 82.496267
## iter 10 value 40.365053
## iter 20 value 19.828875
## iter 30 value 19.532020
## iter 40 value 18.475344
## iter 50 value 17.734886
## final value 17.732462
## convergedY<-predict(iris.nnet,test.data,type="class")
table(test.data[,5],Y)## Y
## C S V
## C 23 0 2
## S 0 25 0
## V 0 0 25主成分分析と重回帰分析を組み合わせる
- 放送大学「心理統計法」第14章 多変量解析のすすめ
- https://statrstart.github.io/2014/09/19/主成分分析/
- https://www.karada-good.net/analyticsr/r-522
library(knitr)
mizumizusisa<-c(3,6,6,6,7,6,5,6,5,5,5,7,7,4,5,5,5,4,6,6,6,7,4,5,5,5,5,6)
namerakasa<-c(6,5,7,1,6,2,5,4,5,6,6,3,3,4,5,3,4,4,5,5,3,5,5,4,5,5,4,3)
simikominohayasa<-c(2,5,5,2,5,6,6,6,5,4,6,6,5,5,6,5,7,5,3,3,2,5,3,4,6,5,4,4)
sarasarakan<-c(2,3,6,7,3,6,3,3,5,5,3,3,7,5,3,5,6,4,3,4,5,1,3,4,4,4,4,7)
sitorikan<-c(5,3,3,1,3,1,5,2,3,5,6,4,1,5,6,3,2,5,5,4,6,7,4,5,6,3,4,1)
turuturukan<-c(5,5,5,2,5,4,4,4,5,4,6,4,6,4,4,2,4,4,5,4,5,4,5,4,4,4,5,3)
yawarakasa<-c(5,6,4,2,4,2,5,1,5,6,7,4,2,5,5,3,4,4,4,5,5,7,5,4,5,3,5,2)
kaorinokonomi<-c(5,6,1,1,1,1,1,1,4,4,1,2,3,1,2,2,2,1,3,5,4,4,3,4,3,5,4,1)
kokotiyosa<-c(6,6,3,2,2,2,3,3,5,5,6,2,3,4,4,4,4,3,4,4,5,4,4,5,4,5,4,1)
d<-data.frame(mizumizusisa,namerakasa,simikominohayasa,sarasarakan,sitorikan,turuturukan,yawarakasa,kaorinokonomi,kokotiyosa,kokotiyosa)
result<- princomp(~mizumizusisa+namerakasa+simikominohayasa+sarasarakan+sitorikan+turuturukan+yawarakasa+kaorinokonomi+kokotiyosa+kokotiyosa)
result## Call:
## princomp(formula = ~mizumizusisa + namerakasa + simikominohayasa +
## sarasarakan + sitorikan + turuturukan + yawarakasa + kaorinokonomi +
## kokotiyosa + kokotiyosa)
##
## Standard deviations:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
## 2.8692000 1.7002112 1.2604634 1.1201672 0.9922762 0.9376346 0.6913266 0.5626037
## Comp.9
## 0.4547934
##
## 9 variables and 28 observations.d<-data.frame(みずみずしさ,なめらかさ,しみこみの速さ,さらさら感“,”しっとり感“,”つるつる感“,”やわらかさ“,”香りの好み“,”心地よさ")
result<- princomp(~“みずみずしさ”+“なめらかさ”+“しみこみの速さ”+“さらさら感”+“しっとり感”+“つるつる感”+“やわらかさ”+“香りの好み”, cor=TRUE, data=d)
- 固有値
result$sdev^2 # component variances## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
## 8.2323089 2.8907181 1.5887680 1.2547746 0.9846121 0.8791586 0.4779324 0.3165229
## Comp.9
## 0.2068370- 主成分負荷量 （テキストとは符号が逆の可能性あり）
library(knitr)
kable(t(result$sdev * t(result$loadings ) )[, 1:3,drop = FALSE],digits=5)| Comp.1 | Comp.2 | Comp.3 | |
|---|---|---|---|
| mizumizusisa | 0.34311 | 0.07882 | 0.17537 |
| namerakasa | -0.92621 | 0.30269 | 0.51943 |
| simikominohayasa | 0.17189 | 0.98593 | 0.68138 |
| sarasarakan | 1.15791 | -0.42973 | 0.12228 |
| sitorikan | -1.46411 | 0.43370 | -0.63268 |
| turuturukan | -0.41739 | 0.00556 | 0.45263 |
| yawarakasa | -1.31321 | 0.17265 | -0.02114 |
| kaorinokonomi | -0.90653 | -1.10821 | 0.35518 |
| kokotiyosa | -1.01107 | -0.43601 | 0.27788 |
- 主成分得点（標準化したもの_テキストとは符号が逆の可能性あり）
result$scores2<-scale(result$scores)
kable(result$scores2[,1:3],digits=4,row.names=T)| Comp.1 | Comp.2 | Comp.3 | |
|---|---|---|---|
| 1 | -1.5034 | -1.4446 | -0.4007 |
| 2 | -0.9860 | -1.2207 | 1.9310 |
| 3 | 0.4232 | 0.7563 | 1.2130 |
| 4 | 2.0952 | -1.2809 | -2.1039 |
| 5 | 0.2808 | 1.2663 | 0.6019 |
| 6 | 1.8291 | 0.3114 | 0.3856 |
| 7 | -0.2468 | 1.6482 | -0.4180 |
| 8 | 1.0551 | 0.8955 | 0.5948 |
| 9 | -0.2572 | -0.6971 | 1.3760 |
| 10 | -0.8443 | -0.5778 | 0.2010 |
| 11 | -1.3066 | 1.5750 | 0.5607 |
| 12 | 0.3998 | 1.0617 | -0.3914 |
| 13 | 1.4407 | -0.9372 | 1.6398 |
| 14 | -0.0421 | 0.7436 | -0.9456 |
| 15 | -0.6502 | 1.2709 | -0.4177 |
| 16 | 0.7633 | -0.1247 | -0.6901 |
| 17 | 0.7504 | 0.4171 | 1.4864 |
| 18 | 0.0970 | 0.9790 | -1.1799 |
| 19 | -0.4974 | -0.2876 | -0.6694 |
| 20 | -0.5077 | -1.2771 | -0.0565 |
| 21 | -0.5807 | -1.4388 | -1.5943 |
| 22 | -1.5693 | 0.7932 | -0.7514 |
| 23 | -0.5612 | -0.4299 | -0.5082 |
| 24 | -0.4482 | -0.7547 | -0.4906 |
| 25 | -0.6202 | 0.7485 | -0.1226 |
| 26 | -0.1403 | -1.0468 | 1.2663 |
| 27 | -0.3593 | -0.6934 | -0.0046 |
| 28 | 1.9861 | -0.2554 | -0.5115 |
- 主成分得点（標準化したもの）上位３つをデータセットdに付け加える
d$PC1<-result$scores2[,1]
d$PC2<-result$scores2[,2]
d$PC3<-result$scores2[,3]
d## mizumizusisa namerakasa simikominohayasa sarasarakan sitorikan turuturukan
## 1 3 6 2 2 5 5
## 2 6 5 5 3 3 5
## 3 6 7 5 6 3 5
## 4 6 1 2 7 1 2
## 5 7 6 5 3 3 5
## 6 6 2 6 6 1 4
## 7 5 5 6 3 5 4
## 8 6 4 6 3 2 4
## 9 5 5 5 5 3 5
## 10 5 6 4 5 5 4
## 11 5 6 6 3 6 6
## 12 7 3 6 3 4 4
## 13 7 3 5 7 1 6
## 14 4 4 5 5 5 4
## 15 5 5 6 3 6 4
## 16 5 3 5 5 3 2
## 17 5 4 7 6 2 4
## 18 4 4 5 4 5 4
## 19 6 5 3 3 5 5
## 20 6 5 3 4 4 4
## 21 6 3 2 5 6 5
## 22 7 5 5 1 7 4
## 23 4 5 3 3 4 5
## 24 5 4 4 4 5 4
## 25 5 5 6 4 6 4
## 26 5 5 5 4 3 4
## 27 5 4 4 4 4 5
## 28 6 3 4 7 1 3
## yawarakasa kaorinokonomi kokotiyosa kokotiyosa.1 PC1 PC2
## 1 5 5 6 6 -1.50339154 -1.4445685
## 2 6 6 6 6 -0.98598687 -1.2207479
## 3 4 1 3 3 0.42318398 0.7563005
## 4 2 1 2 2 2.09522914 -1.2808987
## 5 4 1 2 2 0.28083844 1.2663003
## 6 2 1 2 2 1.82906606 0.3113673
## 7 5 1 3 3 -0.24678340 1.6481514
## 8 1 1 3 3 1.05513903 0.8955204
## 9 5 4 5 5 -0.25715714 -0.6970994
## 10 6 4 5 5 -0.84428810 -0.5777785
## 11 7 1 6 6 -1.30658822 1.5750396
## 12 4 2 2 2 0.39979601 1.0617276
## 13 2 3 3 3 1.44068194 -0.9371894
## 14 5 1 4 4 -0.04209779 0.7435585
## 15 5 2 4 4 -0.65016616 1.2709072
## 16 3 2 4 4 0.76333092 -0.1246799
## 17 4 2 4 4 0.75040020 0.4171061
## 18 4 1 3 3 0.09703209 0.9790010
## 19 4 3 4 4 -0.49738001 -0.2876316
## 20 5 5 4 4 -0.50774166 -1.2770948
## 21 5 4 5 5 -0.58070785 -1.4387566
## 22 7 4 4 4 -1.56925560 0.7932019
## 23 5 3 4 4 -0.56123711 -0.4298600
## 24 4 4 5 5 -0.44815449 -0.7547486
## 25 5 3 4 4 -0.62018010 0.7484690
## 26 3 5 5 5 -0.14033354 -1.0467657
## 27 5 4 4 4 -0.35933807 -0.6934267
## 28 2 1 1 1 1.98608985 -0.2554047
## PC3
## 1 -0.400743191
## 2 1.930953219
## 3 1.213035553
## 4 -2.103916087
## 5 0.601886912
## 6 0.385648535
## 7 -0.417966765
## 8 0.594778196
## 9 1.375981806
## 10 0.200970456
## 11 0.560678993
## 12 -0.391392101
## 13 1.639848549
## 14 -0.945640481
## 15 -0.417733365
## 16 -0.690064294
## 17 1.486350643
## 18 -1.179907029
## 19 -0.669377735
## 20 -0.056524629
## 21 -1.594291858
## 22 -0.751377002
## 23 -0.508185451
## 24 -0.490574893
## 25 -0.122626801
## 26 1.266301274
## 27 -0.004589468
## 28 -0.511522987- AICを用いて変数選択 Model <- lm(心地よさ~PC1+PC2+PC3, data=d) select<-step(Model) library(xtable) print(xtable(select),type=“html”)
Model <- lm(kokotiyosa~PC1+PC2+PC3, data=d)
select<-step(Model)## Start: AIC=-20.83
## kokotiyosa ~ PC1 + PC2 + PC3
##
## Df Sum of Sq RSS AIC
## <none> 9.999 -20.834
## - PC3 1 2.1621 12.161 -17.352
## - PC2 1 5.3229 15.321 -10.883
## - PC1 1 28.6236 38.622 15.005library(xtable)
print(xtable(select),type="html")## <!-- html table generated in R 3.6.1 by xtable 1.8-4 package -->
## <!-- Fri Mar 06 17:19:49 2020 -->
## <table border=1>
## <tr> <th> </th> <th> Estimate </th> <th> Std. Error </th> <th> t value </th> <th> Pr(>|t|) </th> </tr>
## <tr> <td align="right"> (Intercept) </td> <td align="right"> 3.8214 </td> <td align="right"> 0.1220 </td> <td align="right"> 31.33 </td> <td align="right"> 0.0000 </td> </tr>
## <tr> <td align="right"> PC1 </td> <td align="right"> -1.0296 </td> <td align="right"> 0.1242 </td> <td align="right"> -8.29 </td> <td align="right"> 0.0000 </td> </tr>
## <tr> <td align="right"> PC2 </td> <td align="right"> -0.4440 </td> <td align="right"> 0.1242 </td> <td align="right"> -3.57 </td> <td align="right"> 0.0015 </td> </tr>
## <tr> <td align="right"> PC3 </td> <td align="right"> 0.2830 </td> <td align="right"> 0.1242 </td> <td align="right"> 2.28 </td> <td align="right"> 0.0319 </td> </tr>
## </table>- 横浜市立大学方式（practice2)でhtml化
library(stargazer)##
## Please cite as:## Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.2. https://CRAN.R-project.org/package=stargazery<-xtable(select)
stargazer(as.data.frame(y),
type = "html")| Statistic | N | Mean | St. Dev. | Min | Pctl(25) | Pctl(75) | Max |
| Estimate | 4 | 0.658 | 2.176 | -1.030 | -0.590 | 1.168 | 3.821 |
| Std. Error | 4 | 0.124 | 0.001 | 0.122 | 0.124 | 0.124 | 0.124 |
| t value | 4 | 5.436 | 17.795 | -8.289 | -4.753 | 9.541 | 31.329 |
| Pr(> | t| ) | 4 | 0.008 | 0.016 | 0 | 0 | 0.01 | 0 |
変数選択 Start: AIC=-4.37 心地よさ ~ PC1 + PC2 + PC3
Df Sum of Sq RSS AIC
PC3 1 0.0009 18.002 -6.3689 18.001 -4.3703 PC2 1 3.6362 21.637 -1.2186 PC1 1 24.4694 42.470 17.6647 Step: AIC=-6.37 心地よさ ~ PC1 + PC2
Df Sum of Sq RSS AIC
18.002 -6.3689
PC2 1 3.6362 21.638 -3.2175 PC1 1 24.4694 42.471 15.6653 結果
Estimate Std. Error t value Pr(>|t|) (Intercept) 3.8214 0.1604 23.83 0.0000 PC1 -0.9520 0.1633 -5.83 0.0000 PC2 -0.3670 0.1633 -2.25 0.0337
library(ggplot2)
library(GGally)## Warning: package 'GGally' was built under R version 3.6.2## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2##
## Attaching package: 'GGally'## The following object is masked from 'package:dplyr':
##
## nasalibrary(tabplot)## Warning: package 'tabplot' was built under R version 3.6.2## Loading required package: bit## Warning: package 'bit' was built under R version 3.6.2## Attaching package bit## package:bit (c) 2008-2012 Jens Oehlschlaegel (GPL-2)## creators: bit bitwhich## coercion: as.logical as.integer as.bit as.bitwhich which## operator: ! & | xor != ==## querying: print length any all min max range sum summary## bit access: length<- [ [<- [[ [[<-## for more help type ?bit##
## Attaching package: 'bit'## The following object is masked from 'package:data.table':
##
## setattr## The following object is masked from 'package:base':
##
## xor## Loading required package: ff## Warning: package 'ff' was built under R version 3.6.2## Attaching package ff## - getOption("fftempdir")=="C:/Users/721540/AppData/Local/Temp/RtmpsXeoys"## - getOption("ffextension")=="ff"## - getOption("ffdrop")==TRUE## - getOption("fffinonexit")==TRUE## - getOption("ffpagesize")==65536## - getOption("ffcaching")=="mmnoflush" -- consider "ffeachflush" if your system stalls on large writes## - getOption("ffbatchbytes")==342894837.76 -- consider a different value for tuning your system## - getOption("ffmaxbytes")==17144741888 -- consider a different value for tuning your system##
## Attaching package: 'ff'## The following objects are masked from 'package:bit':
##
## clone, clone.default, clone.list## The following objects are masked from 'package:utils':
##
## write.csv, write.csv2## The following objects are masked from 'package:base':
##
## is.factor, is.ordered## Loading required package: ffbase## Warning: package 'ffbase' was built under R version 3.6.2## Registered S3 methods overwritten by 'ffbase':
## method from
## [.ff ff
## [.ffdf ff
## [<-.ff ff
## [<-.ffdf ff##
## Attaching package: 'ffbase'## The following objects are masked from 'package:ff':
##
## [.ff, [.ffdf, [<-.ff, [<-.ffdf## The following objects are masked from 'package:base':
##
## %in%, tabletableplot(d)## Warning in tableplot_checkBins(nBins, max(N, 2)): Setting nBins (100) to number
## of rows (28)
library(psych)## Warning: package 'psych' was built under R version 3.6.2##
## Attaching package: 'psych'## The following object is masked from 'package:bit':
##
## keysort## The following object is masked from 'package:mclust':
##
## sim## The following object is masked from 'package:randomForest':
##
## outlier## The following object is masked from 'package:car':
##
## logit## The following objects are masked from 'package:ggplot2':
##
## %+%, alphapsych::pairs.panels(d)
#ggpairs(d,aes_string(colour="kokotiyosa", alpha=0.5))