※上記の広告は60日以上更新のないWIKIに表示されています。更新することで広告が下部へ移動します。

Conditional Weighted Residual

条件付 (conditional) WRES,すなわち,FOCE の近似に対応した WRES,ということで CWRES と名づけられた.


以下に R のソースを貼っておく.
使い方の詳細な説明はまた後ほど,というか on request で.

S-Plus で動かす場合は,subset の扱いに関して修正が必要.また,sqrtm.cpk() 内の 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/")