# ----------------------------------------------------- rslm()の定義 # Random Sampling Linear Model # dat: 最終列に応答変数をおいたデータ行列 rslm <- function(dat){ nc <- ncol(dat) # 列数の取得 # 変数名が付いていないときには仮の名前を付ける if(is.null(colnames(dat))){ colnames(dat) <- paste("Var", 1:nc, sep="") } dat <- data.frame(dat) # データフレームの形式にする nv <- ncol(dat)-1 # 最終列を除いた列数を取得 v.name <- names(dat) # 全ての変数名の取得 resp.name <- v.name[nv+1] # 応答変数(最終列)名の取得 expl.name <- v.name[1:nv] # 説明変数名の取得 x <- sample(1:nv, sample(1:nv, 1)) # 解析に用いる説明変数をランダムに抽出し、 ana.expl.name <- expl.name[x] # その変数名をana.expl.nameに代入する model <- reformulate(ana.expl.name, resp.name) # 解析するモデル式をmodelへ代入 lm.obj <- lm(model, dat) # lm()オブジェクトをlm.objへ代入 ars <- summary(lm(model, dat))$adj.r.squared # 決定係数を抽出してarsへ代入 # MODEL: モデル式, RESULT: lm()オブジェクト, ARS: 決定係数をリストで返す list(MODEL=model, RESULT=lm.obj, ARS=ars) } # ----------------------------------------------------- rslm.simulation()の定義 # RSLMのシミュレーションを行うための関数 # dat: 最終列に応答変数をおいたデータ行列 # rep.num: rslm()を繰り返す数(rep.num=100くらいで大丈夫) rslm.simulate <- function(dat, rep.num){ ars.dat <- numeric(1) # 決定係数を格納しておく変数 for(i in 1:rep.num){ # 別の関数rslm()を先に定義しておくこと!! out <- rslm(dat) if(rep.num == 1){ # 1回目の繰り返しならば: ars.dat <- out$ARS # ars.datに決定係数を代入する lm.object <- out$RESULT # lm.objectにlm()オブジェクトを代入する cat("model.formula: \n"); print(out$MODEL) # モデル式の表示 cat("Adj.R.Squared: \n"); print(ars.dat) # 決定係数の表示 }else{ # 2回目以降の繰り返しならば: if(ars.dat < out$ARS){ # 前回のR^2よりも大きければ ars.dat <- out$ARS # ars.datに新たなR^2を代入し、 lm.object <- out$RESULT # lm()オブジェクトも更新する cat("model.formula: \n"); print(out$MODEL) # モデル式の表示 cat("Adj.R.Squared: \n"); print(ars.dat) # 決定係数の表示 }else{ # 前回のR^2以下なら next # 次のループへ移動 } # 2つ目のif文の終わり } # 1つ目のif文の終わり } # for文の終わり cat("Complete!! \n") # 終了の合図を表示 return(lm.object) # 結局、R^2が最大であったlm()オブジェクトを返す } # ----------------------------------------------------- atehame()の定義 # 欠損値データに対しての当てはめ値を返す関数 # rv: 応答変数(response variables) # prd: vglm()によって得られた予測確率 atehame <- function(rv, prd){ size <- nrow(prd) # ad.vglm.prdの列数 atv <- numeric(size) # 当てはめ値を代入していく変数 for(i in 1:size){ # 列数=sizeの数だけ繰り返す sr <- prd[i, ] == max(prd[i, ]) # atv[i] <- (1:length(sr))[sr == TRUE] # } tbl <- table(abs(rv - atv)) # 判別的中数・誤判別数 prob <- tbl / sum(tbl) # 判別的中率・誤判別率 list(FITTED.VALUE = atv, TABLE = tbl, PROBABILITY = prob) } # ----------------------------------------------------- mlri2()の定義 −本番用− # Multinomial Logistic Regression Imputation # ***注意*** あらかじめatehame()を定義しておく必要がある # dat: 最終列に応答変数を置いたデータ行列 # sn: 欠損値データとみなして抽出する数(データ行数) # rslm.num: rslm.simulate()での繰り返し数 mlri2 <- function(dat, rslm.num){ colnum <- ncol(dat) # 引数datの列数を取得 colnames(dat) <- paste("Var", 1:colnum, sep="") # 列ラベルをつける dat <- data.frame(dat) # データフレームに変換 # 欠損値がある行と欠損値がない行とを分離する mis.rows <- complete.cases(dat) # 欠損値のある行がFALSEとされる ana.dat <- dat[mis.rows, ] # TRUEである行のみを抽出する: 解析用データ no.mis.rows <- mis.rows==FALSE # 欠損値のない行がFALSEとされる mis.dat <- dat[no.mis.rows, ] # TRUEである行のみを抽出する: 欠損値データ # rslm.simulate()の出力が膨大になるのでテキストファイルに出力する sink("rslm.out.txt") rslm.result <- rslm.simulate(ana.dat, rslm.num) # 引数rslm.numで指定された数を繰り返す sink() cat("rslm.out.txtとしてワークディレクトリに出力されました \n") ad.lm.prd <- predict(rslm.result) # 解析用データの重回帰予測値 ad.Y <- ana.dat[, colnum] # 解析用データの応答変数 # もしVGAMパッケージが読み込まれていなければVGAMを読み込む if(any(search() == "package:VGAM") == FALSE) library(VGAM) vglm.result <- vglm(as.factor(ad.Y) ~ ad.lm.prd, multinomial) # 多項ロジットモデル ad.vglm.prd <- predict(vglm.result, type="response") # 解析用データの予測確率 ad.res <- atehame(ad.Y, ad.vglm.prd) # 解析用データへの当てはめ値 md.lm.prd <- predict(rslm.result, mis.dat) # 欠損値データの重回帰予測値 md.vglm.prd <- predict(vglm.result, # 欠損値データの予測確率 list(ad.lm.prd=md.lm.prd), type="response") md.res <-atehame(mis.dat[,colnum], md.vglm.prd) # 欠損値データへの当てはめ値 list(MULTIPLE = rslm.result, LOGIT = vglm.result, ANA.MR.PREDICT = ad.lm.prd, ANA.LR.PREDICT = ad.vglm.prd, MIS.MR.PREDICT = md.lm.prd, MIS.LR.PREDICT = md.vglm.prd, ANA.RESULT = ad.res, MIS.RESULT = md.res) }