×
[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。
オブジェクトにコメントをつける。データフレームや解析結果につけるようにしよう
df <- iris
comment(df) <- c("iris dataframe", "日本語はどうだろう")
comment(df)
[1] "iris data frame" "日本語はどうだろう"
# 日本語もいけるようだ。ちなみにdfだけだと何もでない
作業スペースの保存
save.image("working.RData")
# コンソールの保存はファイルの保存を選択しないとできないみたい
出力をファイルに記録する
sink("log.txt", append=TRUE, split=TRUE)
# 出力を作業ディレクトリのlog.txtというテキストファイルに保存する
# append=TRUEは追加モード (上書きしない) 。split=TRUEはコンソールにも出力する。
データフレームをテキストで保存する
write.table(df, file="iris.txt", sep="\t", row.names=FALSE, quote=FALSE)
# sep="\t"はタブ区切り、row.names=FALSEは行番号を削除、quote=FALSEはダブルクォーテーション削除
テキストファイルからデータフレームを読み込む
df <- read.table("Data/iris.txt", header=T)
？？
# sink関数でtype = c("output", "message")という引数も指定できるが何かわからない。
# 公式だとtypeは以下のようなことらしい。stdoutとstderrって何だ？UNIXの何か？
Normal R output (to connection stdout)) is diverted by the default type = "output". Only prompts and (most) messages continue to appear on the console. Messages sent to stderr() (including those from message, warning and stop) can be diverted by sink(type = "message")
df <- iris
comment(df) <- c("iris dataframe", "日本語はどうだろう")
comment(df)
[1] "iris data frame" "日本語はどうだろう"
# 日本語もいけるようだ。ちなみにdfだけだと何もでない
作業スペースの保存
save.image("working.RData")
# コンソールの保存はファイルの保存を選択しないとできないみたい
出力をファイルに記録する
sink("log.txt", append=TRUE, split=TRUE)
# 出力を作業ディレクトリのlog.txtというテキストファイルに保存する
# append=TRUEは追加モード (上書きしない) 。split=TRUEはコンソールにも出力する。
データフレームをテキストで保存する
write.table(df, file="iris.txt", sep="\t", row.names=FALSE, quote=FALSE)
# sep="\t"はタブ区切り、row.names=FALSEは行番号を削除、quote=FALSEはダブルクォーテーション削除
テキストファイルからデータフレームを読み込む
df <- read.table("Data/iris.txt", header=T)
？？
# sink関数でtype = c("output", "message")という引数も指定できるが何かわからない。
# 公式だとtypeは以下のようなことらしい。stdoutとstderrって何だ？UNIXの何か？
Normal R output (to connection stdout)) is diverted by the default type = "output". Only prompts and (most) messages continue to appear on the console. Messages sent to stderr() (including those from message, warning and stop) can be diverted by sink(type = "message")
PR
順序ロジスティック
従属変数が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)))
先に結論: 単にクロス集計表をつくるなら、Excelのピボットテーブルの方が楽
table, xtabs, ftableがある
組み込みデータフレームesophを使う
"agegp","alcgp","tobgp","ncases","ncontrols" の5列
table(esoph[c(1,2,3)])
1(agegp) が縦軸、2 (alcgp) を横軸にして、3 (tobgp) の群の数だけクロス集計表ができる
2元の集計表が複数ほしいときに使おう
xtabs(~ agegp + alcgp + tobgp, esoph)
上のtableと同じ表を返す
チルダの前にnumericの列をおくと、グループごとに数値を合計した結果が返される
xtabs(ncases ~ agegp + alcgp, esoph)
ftable(esoph[,1], esoph[,2], esoph[,3])
3元 (以上) の集計表をつくりたいときはこれ
一番最後の変数 (esoph[,3]) が横軸 (上) にきて、あとは左から並ぶ
row.vars=3, row.vars=c(1,3)...とかで左 (同様にcol.var...) にくる変数を選べる
小計欄が必要なときは
## 先にxtabsで2元のグラフをつくり、
xtbs <- xtabs(~ agegp + alcgp + tobgp, esoph)
## 合計欄を追加して
xmgn <- addmargins(xtbs)
## ftableに渡す
ftable(xmgn)
## ftableにaddmargins関数を使うとなぜか変数名が消えて行列になるのでやめよう
明らかにExcelのピボットテーブルの方が楽
ただ、xtabs, table関数の結果は"table"というクラスになり (xtabsは"xtabs"というのももつ)
summary(xtbs) としたり、chisq.test(xtbs) でカイ二乗検定を行えるのでその点は便利
chisq.testは2次元表しか対応してないため、3次元表も2次元行列と解釈してしまうが、summaryは3次元表としてカイ二乗検定をしてくれるそうだ (xtabsヘルプより)
(もちろんこのとき合計欄をつけてはいけない)
青木先生の関数を使うとファイ係数、一致係数、クラメール係数などが算出できる。感謝
phi, contingency, cramer, chisq
vcdパッケージについては今度調べよう
中澤先生のサイトに使用例あり。感謝
table, xtabs, ftableがある
組み込みデータフレームesophを使う
"agegp","alcgp","tobgp","ncases","ncontrols" の5列
table(esoph[c(1,2,3)])
1(agegp) が縦軸、2 (alcgp) を横軸にして、3 (tobgp) の群の数だけクロス集計表ができる
2元の集計表が複数ほしいときに使おう
xtabs(~ agegp + alcgp + tobgp, esoph)
上のtableと同じ表を返す
チルダの前にnumericの列をおくと、グループごとに数値を合計した結果が返される
xtabs(ncases ~ agegp + alcgp, esoph)
ftable(esoph[,1], esoph[,2], esoph[,3])
3元 (以上) の集計表をつくりたいときはこれ
一番最後の変数 (esoph[,3]) が横軸 (上) にきて、あとは左から並ぶ
row.vars=3, row.vars=c(1,3)...とかで左 (同様にcol.var...) にくる変数を選べる
小計欄が必要なときは
## 先にxtabsで2元のグラフをつくり、
xtbs <- xtabs(~ agegp + alcgp + tobgp, esoph)
## 合計欄を追加して
xmgn <- addmargins(xtbs)
## ftableに渡す
ftable(xmgn)
## ftableにaddmargins関数を使うとなぜか変数名が消えて行列になるのでやめよう
明らかにExcelのピボットテーブルの方が楽
ただ、xtabs, table関数の結果は"table"というクラスになり (xtabsは"xtabs"というのももつ)
summary(xtbs) としたり、chisq.test(xtbs) でカイ二乗検定を行えるのでその点は便利
chisq.testは2次元表しか対応してないため、3次元表も2次元行列と解釈してしまうが、summaryは3次元表としてカイ二乗検定をしてくれるそうだ (xtabsヘルプより)
(もちろんこのとき合計欄をつけてはいけない)
青木先生の関数を使うとファイ係数、一致係数、クラメール係数などが算出できる。感謝
phi, contingency, cramer, chisq
vcdパッケージについては今度調べよう
中澤先生のサイトに使用例あり。感謝
データフレームの列 (変数ごと) に記述統計を出したい
最小値，第1 四分位点，中央値，平均，第3 四分位点，最大値を求める
summary(dat) # NAは自動的に除かれる
不偏標準偏差、不偏分散、標準誤差、データ数 (NAこみ、NAなし、NA) 、歪度、尖度、標本標準偏差、標本分散を求める
##メモ: na.omit(x)はNAを除いてベクトルの要素をそのまま返す。!is.na(x)はTRUE/FALSEを返す。だから数を数えたいときは!is.naを使う。
sd(dat, na.rm=T) # 不偏標準偏差
sd(dat, na.rm=T)^2 # 不偏分散
sapply(dat, function(x) sqrt(var(as.vector(na.omit(x)))/length(na.omit(x)))) 標準誤差
sapply(dat, length) # データ数、NAこみ
sapply(dat, function(x) length(na.omit(x))) # データ数、NAなし
sapply(dat, function(x) sum(is.na(x))) # NAの数
sapply(dat, function(x) mean((na.omit(x)-mean(na.omit(x)))^3)/(sd(na.omit(x))^3)) # 歪度。データにNAがあった場合、省いて計算したいできた
sapply(dat, function(x) mean((na.omit(x)-mean(na.omit(x)))^4)/(sd(na.omit(x))^4)) # 尖度。上に同じ
sapply(dat, function(x) sqrt(var(na.omit(x))*(length(na.omit(x))-1)/length(na.omit(x)))) # 標本標準偏差
sapply(dat, function(x) var(na.omit(x))*(length(na.omit(x))-1)/length(na.omit(x))) # 標本分散。#なぜかわからんがvarは因子変数にも適用できるようだ
# var(x, na.rm=T) でも同じ
グループ別の統計量
aggregate(iris[1:4], iris[5], mean, na.rm = T)
その他
rowMeans(dat, na.rm=TRUE) # 行ごとの平均、NAを抜く
rowSums(dat, na.rm=TRUE) # 行ごとの合計、NAを抜く
apply(dat, 1, function(x) sum(is.na(x))) #行ごとのNA数
sum(complete.cases(dat)) # データフレーム全体で欠損値の無いケース数を調べる
sum((na.omit(x)-mean(na.omit(x)))^2) # 平方和
最小値，第1 四分位点，中央値，平均，第3 四分位点，最大値を求める
summary(dat) # NAは自動的に除かれる
不偏標準偏差、不偏分散、標準誤差、データ数 (NAこみ、NAなし、NA) 、歪度、尖度、標本標準偏差、標本分散を求める
##メモ: na.omit(x)はNAを除いてベクトルの要素をそのまま返す。!is.na(x)はTRUE/FALSEを返す。だから数を数えたいときは!is.naを使う。
sd(dat, na.rm=T) # 不偏標準偏差
sd(dat, na.rm=T)^2 # 不偏分散
sapply(dat, function(x) sqrt(var(as.vector(na.omit(x)))/length(na.omit(x)))) 標準誤差
sapply(dat, length) # データ数、NAこみ
sapply(dat, function(x) length(na.omit(x))) # データ数、NAなし
sapply(dat, function(x) sum(is.na(x))) # NAの数
sapply(dat, function(x) mean((na.omit(x)-mean(na.omit(x)))^3)/(sd(na.omit(x))^3)) # 歪度。データにNAがあった場合、省いて計算
sapply(dat, function(x) mean((na.omit(x)-mean(na.omit(x)))^4)/(sd(na.omit(x))^4)) # 尖度。上に同じ
sapply(dat, function(x) sqrt(var(na.omit(x))*(length(na.omit(x))-1)/length(na.omit(x)))) # 標本標準偏差
sapply(dat, function(x) var(na.omit(x))*(length(na.omit(x))-1)/length(na.omit(x))) # 標本分散。#なぜかわからんがvarは因子変数にも適用できるようだ
# var(x, na.rm=T) でも同じ
グループ別の統計量
aggregate(iris[1:4], iris[5], mean, na.rm = T)
その他
rowMeans(dat, na.rm=TRUE) # 行ごとの平均、NAを抜く
rowSums(dat, na.rm=TRUE) # 行ごとの合計、NAを抜く
apply(dat, 1, function(x) sum(is.na(x))) #行ごとのNA数
sum(complete.cases(dat)) # データフレーム全体で欠損値の無いケース数を調べる
sum((na.omit(x)-mean(na.omit(x)))^2) # 平方和
