×
[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。
## ロジスティック回帰の一般的注意
独立変数に名義尺度を使った場合、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)))
PR
Comment
Trackback
Trackback URL
Comment form
