「CWRES」の編集履歴（バックアップ）一覧に戻る
CWRES - (2006/06/26 (月) 22:40:36) のソース
*Conditional Weighted Residual
条件付 (conditional) WRES，すなわち，FOCE の近似に対応した WRES，ということで CWRES と名づけられた．
-[[Xpose のサイト>http://xpose.sourceforge.net/]]
-[[今年 (2006) の PAGE での発表>http://www.page-meeting.org/?abstract=1001]]
-[[関連記事>http://blog.goo.ne.jp/hkasai/e/ab69121c9319ece3655d3fe9b1215b2d]] (コントロールファイルの記述方法，等)
以下に R のソースを貼っておく．
使い方の詳細な説明はまた後ほど，というか on request で．
S-Plus で動かす場合は，subset の扱いに関して修正が必要．また，sqrtmcpk() 内の eigen(x) を eigen(x, symmetric=T) としてください．
----
### Conditional WRES (CWRES) の計算 ###
compute.cwres.cpk <-
function(
tab.cwres.fname, # G, HH の情報の入った TABLE ファイル
# このファイルには以下の列が必要
# ID, MDV, DV, IPRED, G11, G21, ..., H11, H21, ...
# (MDV がない場合でも計算可能だが面倒なので未対応)
n.theta, # THETA の数
n.eta, # ETA の数
n.eps=1, # EPS の数
ID="ID", # 被験者番号の列名
DV="DV", # DV の列名
IPRE="IPRE", # IPRED の列名
tab.par.fname="par.txt", # INFN 出力のパラメータ推定値ファイル
tab.eta.fname="eta.txt", # INFN 出力の ETA POSTHOC 推定値ファイル
folder # 実行フォルダ名 (上記ファイルが入っているフォルダ)
################################################################################
### ######
### 結果: 元の TABLE ファイルの MDV=0 のレコードについて計算した CWRES を ######
### データセットの右端に付与して返す ######
### ######
################################################################################
) {
### フォルダ名の最後に '/' がなかったら付与する
if (folder[length(folder)] != '/')
folder <- paste(folder, '/', sep="")
tab.cwres.all <-
read.table(paste(folder, tab.cwres.fname, sep=""), header=T, skip=1)
### MDV 列があるなら MDV=0 のレコードのみを選択 ###
if (!is.null(tab.cwres.all$MDV))
tab.cwres <- tab.cwres.all[tab.cwres.all$MDV==0,]
### ないときはエラー。****** 要 REFINE ****** ###
else
stop("MDV 列を TABLE に含めてください．\n")
### MDV=0 のインデックス ###
idx.MDV <- (tab.cwres.all$MDV == 0)
### ID 番号を取り出す ###
ID.uniq <- tab.cwres[!duplicated(tab.cwres[ID]), ID]
tab.par <- scan(paste(folder, tab.par.fname, sep=""))
obj <- tab.par[1] ### 目的関数値 ###
iere <- tab.par[2]; ierc <- tab.par[3] ### EST, COV のエラーコード ###
### THETA の推定値 ###
THETA <- tab.par[seq(4, length=n.theta)]
### OMEGA の推定値 ###
OMEGA <- matrix(tab.par[seq(4+2*n.theta, length=n.eta^2)], ncol=n.eta)
### SIGMA の推定値 ###
SIGMA <- matrix(tab.par[seq(4+2*n.theta+2*n.eta^2, length=n.eps^2)], ncol=n.eps)
### ETA の POSTHOC 推定値 ###
### 被験者数 x ETA の数、の形式 ###
tab.eta <- read.table(paste(folder, tab.eta.fname, sep=""), header=F)
### ETA1, ETA2, ... と変数名をつける
names(tab.eta) <- paste("ETA", seq(ncol(tab.eta)), sep="")
tab.eta <- cbind(ID=ID.uniq, tab.eta)
### 被験者ごとに CWRES を計算してまとめる ###
cwres <- unlist(by(tab.cwres, tab.cwres[ID],
FUN=compute.cwres.indiv,
THETA=THETA, OMEGA=OMEGA, SIGMA=SIGMA, ETA=tab.eta,
ID, DV, IPRE
))
#cbind(tab.cwres, CWRES=cwres)
tab.cwres.all$CWRES <- rep(0.0, nrow(tab.cwres.all))
tab.cwres.all$CWRES[idx.MDV] <- cwres
tab.cwres.all
}
### 被験者ごとの CWRES 計算 ###
compute.cwres.indiv <-
function(data, THETA, OMEGA, SIGMA, ETA, ID, DV, IPRE) {
n.eta <- ncol(OMEGA)
G.colnames <- paste("G", seq(1, n.eta), "1", sep="")
### G11, G21, ... のみを抽出 ###
G.matrix <- as.matrix(subset(data, select=G.colnames))
n.eps <- ncol(SIGMA)
H.colnames <- paste("H", seq(1, n.eps), "1", sep="")
### H11, H21, ... のみを抽出 ###
H.matrix <- as.matrix(subset(data, select=H.colnames))
### 対角項のみのベクトルになる ###
H.sigma.H <- diag(H.matrix %*% SIGMA %*% t(H.matrix))
### 残差の分散共分散行列 ###
COV.indiv <- G.matrix %*% OMEGA %*% t(G.matrix) + diag(H.sigma.H)
### ETA の POSTHOC 推定値 ###
ETA.indiv <- ETA[ETA$ID==data[1, ID],]
CIPRE <- data[IPRE] -
G.matrix %*% t(as.matrix(ETA.indiv[, -1])) ### ETA.indiv の第 1 列は ID ###
CIRES <- as.matrix(data[DV] - CIPRE)
### COV.indiv^(-1/2) %*% CIRES を求める ###
CWRES.indiv <- solve(sqrtm.cpk(COV.indiv), CIRES)
CWRES.indiv
}
### 行列 x の平方根: x = x.half %*% x.half ###
sqrtm.cpk <-
function(x) {
x.eigen <- eigen(x)
H <- x.eigen$vectors
ev <- x.eigen$values
ev[ev < 0] <- 0
Dlambda.half <- diag(sqrt(ev))
x.half <- H %*% Dlambda.half %*% t(H)
x.half
}
#compute.cwres.cpk("tabcwres.txt", n.theta=2, n.eta=2, n.eps=1, ID="SID", folder="d:/nmv/run/cwres/")
