×
[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。
正準判別分析をやってみた2
(正準判別分析をやってみた1)
データはエヴェリット (2007) の頭蓋骨データ
青木先生のcandis関数
rescandis <- candis(skulls[2:5],skulls[1])
# 結果print(rescandis)
基本はhttp://aoki2.si.gunma-u.ac.jp/R/candis.htmlを参照
lda関数、エヴェリット (2007) と関連するところだけ列挙
判別係数 #ldaのCoefficients of linear discriminants、エヴェリット (2007) の正準係数
axis 1 axis 2 axis 3 axis 4
MB -0.126676 0.038738 0.092768 -0.14883986
- 中略 -
NH -0.082851 -0.077293 -0.294589 -0.06685888
constant 2.283887 -22.530493 4.409749 36.07126158
判別結果 #エヴェリット (2007) のプラグイン推定による分類表
classification
group
c1850BC c200BC c3300BC c4000BC cAD150
c1850BC 15 2 4 4 5
- 中略 -
cAD150 4 9 4 2 11
$pooled.cov #エヴェリット (2007) の共分散構造：等分散性
MB BH BL NH
MB 21.110805 0.036782 0.07908 2.0090
BH 0.036782 23.484598 5.20000 2.8451
BL 0.079080 5.200000 24.17908 1.1333
NH 2.008966 2.845057 1.13333 10.1526
その他
ans <- predict(reslda) #lda関数による判別の結果を返す
candis関数では以下に一致
$p.Bayes はベイズ確率。predict(reslda)の$posterior。判別後のケースごとの確率
$can.score は判別得点。predict(reslda)の$xと一致
$classificationは分類。predict(reslda)の$class (だと思う)
# groupは元の群、classificationは判別関数による分類
candis関数のprintについて
printで判別結果がうまくでないときがある。resultオブジェクト生成の問題のようだ
(正準判別分析をやってみた1)
データはエヴェリット (2007) の頭蓋骨データ
青木先生のcandis関数
rescandis <- candis(skulls[2:5],skulls[1])
# 結果print(rescandis)
基本はhttp://aoki2.si.gunma-u.ac.jp/R/candis.htmlを参照
lda関数、エヴェリット (2007) と関連するところだけ列挙
判別係数 #ldaのCoefficients of linear discriminants、エヴェリット (2007) の正準係数
axis 1 axis 2 axis 3 axis 4
MB -0.126676 0.038738 0.092768 -0.14883986
- 中略 -
NH -0.082851 -0.077293 -0.294589 -0.06685888
constant 2.283887 -22.530493 4.409749 36.07126158
判別結果 #エヴェリット (2007) のプラグイン推定による分類表
classification
group
c1850BC c200BC c3300BC c4000BC cAD150
c1850BC 15 2 4 4 5
- 中略 -
cAD150 4 9 4 2 11
$pooled.cov #エヴェリット (2007) の共分散構造：等分散性
MB BH BL NH
MB 21.110805 0.036782 0.07908 2.0090
BH 0.036782 23.484598 5.20000 2.8451
BL 0.079080 5.200000 24.17908 1.1333
NH 2.008966 2.845057 1.13333 10.1526
その他
ans <- predict(reslda) #lda関数による判別の結果を返す
candis関数では以下に一致
$p.Bayes はベイズ確率。predict(reslda)の$posterior。判別後のケースごとの確率
$can.score は判別得点。predict(reslda)の$xと一致
$classificationは分類。predict(reslda)の$class (だと思う)
# groupは元の群、classificationは判別関数による分類
candis関数のprintについて
printで判別結果がうまくでないときがある。resultオブジェクト生成の問題のようだ
PR
正準判別分析をやってみた
データは
B.エヴェリット (2007) のp108-110 表5.8 頭蓋骨データ
リンク中のchap5skulls
従属変数になる群データはEPOCH (1列目) 、これをMB, BH, BL, NH (2-5列目) という連続変量で予測する
MASSパッケージのlda関数
# 線型判別の関数だが、3群以上あると正準判別をやってくれるらしい
reslda <- lda(skulls[,2:5], skulls[,1])
# これだけ？
print(reslda) # 結果
Call: lda(skulls[, 2:5], skulls[, 1])
# 事前確率
Prior probabilities of groups:
c4000BC c3300BC c1850BC c200BC cAD150
0.2 0.2 0.2 0.2 0.2
# 群内の平均値
Group means:
MB BH BL NH
c4000BC 131.3667 133.6000 99.16667 50.53333
-中略-
cAD150 136.1667 130.3333 93.50000 51.36667
# 判別係数 (正準係数とか正準判別係数とも)
Coefficients of linear discriminants:
LD1 LD2 LD3 LD4
MB 0.12667629 0.03873784 0.09276835 0.1488398644
-中略- BH, BLの結果
NH 0.08285128 -0.07729281 -0.29458931 0.0668588797
Proportion of trace:
LD1 LD2 LD3 LD4
0.8823 0.0809 0.0326 0.0042
結果がうまくhtmlででないなぁ…
データは
B.エヴェリット (2007) のp108-110 表5.8 頭蓋骨データ
リンク中のchap5skulls
従属変数になる群データはEPOCH (1列目) 、これをMB, BH, BL, NH (2-5列目) という連続変量で予測する
MASSパッケージのlda関数
# 線型判別の関数だが、3群以上あると正準判別をやってくれるらしい
reslda <- lda(skulls[,2:5], skulls[,1])
# これだけ？
print(reslda) # 結果
Call: lda(skulls[, 2:5], skulls[, 1])
# 事前確率
Prior probabilities of groups:
c4000BC c3300BC c1850BC c200BC cAD150
0.2 0.2 0.2 0.2 0.2
# 群内の平均値
Group means:
MB BH BL NH
c4000BC 131.3667 133.6000 99.16667 50.53333
-中略-
cAD150 136.1667 130.3333 93.50000 51.36667
# 判別係数 (正準係数とか正準判別係数とも)
Coefficients of linear discriminants:
LD1 LD2 LD3 LD4
MB 0.12667629 0.03873784 0.09276835 0.1488398644
-中略- BH, BLの結果
NH 0.08285128 -0.07729281 -0.29458931 0.0668588797
Proportion of trace:
LD1 LD2 LD3 LD4
0.8823 0.0809 0.0326 0.0042
結果がうまくhtmlででないなぁ…
順序ロジスティック
従属変数が3分類以上で順序尺度のとき
MASSパッケージのpolr関数
累積ロジスティック (比例ロジスティック) モデルに従う
多項ロジスティック
従属変数が3分類以上で名義尺度のとき
nnetパッケージのmultinom関数
多重ロジスティック
独立変数が複数あるとき。つまり普通ロジスティックといえばこれ
基本パッケージのglm, nls関数、
Designパッケージのlrm関数などを使う
glmは尤離度 (deviance, 逸脱度) 、nlsは残差を最小化するので結果が違う
lrmはhttp://stat.cmu.edu/S/Harrell/help/Design/html/lrm.html
従属変数が3分類以上で順序尺度のとき
MASSパッケージのpolr関数
累積ロジスティック (比例ロジスティック) モデルに従う
多項ロジスティック
従属変数が3分類以上で名義尺度のとき
nnetパッケージのmultinom関数
多重ロジスティック
独立変数が複数あるとき。つまり普通ロジスティックといえばこれ
基本パッケージのglm, nls関数、
Designパッケージのlrm関数などを使う
glmは尤離度 (deviance, 逸脱度) 、nlsは残差を最小化するので結果が違う
lrmはhttp://stat.cmu.edu/S/Harrell/help/Design/html/lrm.html
## ロジスティック回帰の一般的注意
独立変数に名義尺度を使った場合、Rではfactor(変数, levels=c("0", "1"...)) で変換する。分析時にはこれが自動的にダミー変数化される。なので、カテゴリーをむやみにカテゴリーを増やすとケース数より多くなっちゃうかもしれない
http://aoki2.si.gunma-u.ac.jp/lecture/Regression/mreg/dummy-variable/dummy.html
http://qdai.way-nifty.com/qjes/2008/01/r_85a3.html
http://www.aichi-gakuin.ac.jp/~chino/multivar/chapter4/sec4.html
http://www.ibaraki-kodomo.com/toukei/logis.html
Designパッケージのlrm関数を使うと下みたいにいろいろ自分でいろいろ計算しなくていいので楽…らしい。詳しくはそのうち調べよう
# 大本
result <- glm(y~x1+x2, dat, family=binomial)
#因子変数は順序尺度変数にする。順序尺度の一番下 (levels=c("0", "1", "2")...とすると0のケースはリファレンスカテゴリとなり、0と比べ、2のケース、3のケースは何倍のオッズになるかがわかる。結局は比の問題なの で何にしようがオッズ比のp値は変わらないが
#尤度比検定
#独立変数がないモデルをつくり、比較する。検定しているのは-2*対数尤度
result0 <- glm(y ~ 1, dat, family=binomial)
anova(result0, result, test="Chisq")
#個々の対数尤度はlogLik
logLik(result)
logLik(result0)
#回帰式
summary(result)
##信頼区間と標準化偏回帰係数は自分で計算する
###Std.Error*1.96 のプラスマイナスで信頼区間
result.se <- summary(result)$coefficients[, "Std. Error"]
result.uci <- coef(result)+(result.se*1.96) #上限
result.lci <- coef(result)-(result.se*1.96) #下限
print(result.uci)
print(result.lci)
###標準化偏回帰係数は元の独立変数の標準偏差を回帰係数にかける
# 独立変数の標準偏差を求め、ベクトルにする。
xsd <- c(0, sd(dat$x1), sd(dat$x2))
# 編回帰係数に標準偏差を掛けると標準化編回帰係数に
print(beta <- coef(result)*xsd)
これ嘘。そのうち調べよう
#Wald統計量
## (偏回帰係数/その標準誤差)^2 で算出できる。
## p値はsummaryで出る回帰係数と全く同じ
### この統計量、何の意味が…
result.se <- summary(result)$coefficients[, "Std. Error"]
result.wald <- ((result$coef)/(result.se))^2
summary(result)の Coefficients:のz valueで出てた
#オッズ比と信頼区間
print(exp(coef(result))) #オッズ比
# 京都大学のサイトに便利な関数
LOG <- function(tmp) {
lreg.coeffs <- coef(summary(tmp))
ORlow <- exp(lreg.coeffs[,1]-qnorm(0.975)*lreg.coeffs[,2])
OR <- exp(lreg.coeffs[,1])
ORup <- exp(lreg.coeffs[,1]+qnorm(0.975)*lreg.coeffs[,2])
P <-pnorm(abs(lreg.coeffs[,3]),lower.tail=F)*2
tmp <- cbind(OR, ORlow, ORup, P)
return(tmp)
}
LOG(result)
#理論値 (確率)
result$fitted.values
#ロジット予測式の分母
result$linear.predictors
#Nagelkerke's R二乗
# 中澤先生のサイトに便利な関数
NagelkerkeR2 <- function(rr,n) {(1-exp((rr$dev-rr$null)/n))/(1-exp(-rr$null/n)) }
print(NagelkerkeR2(result,nrow(dat)))
独立変数に名義尺度を使った場合、Rではfactor(変数, levels=c("0", "1"...)) で変換する。分析時にはこれが自動的にダミー変数化される。なので、カテゴリーをむやみにカテゴリーを増やすとケース数より多くなっちゃうかもしれない
http://aoki2.si.gunma-u.ac.jp/lecture/Regression/mreg/dummy-variable/dummy.html
http://qdai.way-nifty.com/qjes/2008/01/r_85a3.html
http://www.aichi-gakuin.ac.jp/~chino/multivar/chapter4/sec4.html
http://www.ibaraki-kodomo.com/toukei/logis.html
Designパッケージのlrm関数を使うと下みたいにいろいろ自分でいろいろ計算しなくていいので楽…らしい。詳しくはそのうち調べよう
# 大本
result <- glm(y~x1+x2, dat, family=binomial)
#因子変数は順序尺度変数にする。順序尺度の一番下 (levels=c("0", "1", "2")...とすると0のケースはリファレンスカテゴリとなり、0と比べ、2のケース、3のケースは何倍のオッズになるかがわかる。結局は比の問題なの で何にしようがオッズ比のp値は変わらないが
#尤度比検定
#独立変数がないモデルをつくり、比較する。検定しているのは-2*対数尤度
result0 <- glm(y ~ 1, dat, family=binomial)
anova(result0, result, test="Chisq")
#個々の対数尤度はlogLik
logLik(result)
logLik(result0)
#回帰式
summary(result)
##信頼区間と標準化偏回帰係数は自分で計算する
###Std.Error*1.96 のプラスマイナスで信頼区間
result.se <- summary(result)$coefficients[, "Std. Error"]
result.uci <- coef(result)+(result.se*1.96) #上限
result.lci <- coef(result)-(result.se*1.96) #下限
print(result.uci)
print(result.lci)
# 独立変数の標準偏差を求め、ベクトルにする。
xsd <- c(0, sd(dat$x1), sd(dat$x2))
# 編回帰係数に標準偏差を掛けると標準化編回帰係数に
print(beta <- coef(result)*xsd)
これ嘘。そのうち調べよう
#Wald統計量
## (偏回帰係数/その標準誤差)^2 で算出できる。
## p値はsummaryで出る回帰係数と全く同じ
### この統計量、何の意味が…
result.wald <- ((result$coef)/(result.se))^2
summary(result)の Coefficients:のz valueで出てた
#オッズ比と信頼区間
print(exp(coef(result))) #オッズ比
# 京都大学のサイトに便利な関数
LOG <- function(tmp) {
lreg.coeffs <- coef(summary(tmp))
ORlow <- exp(lreg.coeffs[,1]-qnorm(0.975)*lreg.coeffs[,2])
OR <- exp(lreg.coeffs[,1])
ORup <- exp(lreg.coeffs[,1]+qnorm(0.975)*lreg.coeffs[,2])
P <-pnorm(abs(lreg.coeffs[,3]),lower.tail=F)*2
tmp <- cbind(OR, ORlow, ORup, P)
return(tmp)
}
LOG(result)
#理論値 (確率)
result$fitted.values
#ロジット予測式の分母
result$linear.predictors
#Nagelkerke's R二乗
# 中澤先生のサイトに便利な関数
NagelkerkeR2 <- function(rr,n) {(1-exp((rr$dev-rr$null)/n))/(1-exp(-rr$null/n)) }
print(NagelkerkeR2(result,nrow(dat)))
