「lmeとnlme」の編集履歴（バックアップ）一覧はこちら

「lmeとnlme」(2008/04/05 (土) 16:45:11) の最新版変更点

追加された行は緑色になります。

削除された行は赤色になります。

*変量効果、個体内誤差のモデリング

 require(nlme)
 
 set.seed(2801)
 
 # lme
 
 n <- 100
 m <- 5
 
 data.lme <- data.frame(
   ID=rep(seq(n), each=m),
   TIME=runif(n*m, max=10),
   k0=rep(10*exp(rnorm(n, sd=0.2)), each=m),
   k1=rep(1*exp(rnorm(n, sd=0.2)), each=m)
 )
 data.lme$pred <- with(data.lme, k0+k1*TIME)
 data.lme$y <- data.lme$pred * exp(rnorm(nrow(data.lme), sd=0.1))
 
 summary(fm1 <- lm(y ~ TIME, data.lme))
 
 # 切片のみ。ID 単位で個体間変動を考える
 summary(fm2 <- lme(y ~ TIME, data.lme, random=list(ID=~1)))
 
 # 切片と傾き、独立
 summary(fm3 <- update(fm2, random=list(ID=pdDiag(~1+TIME))))
 
 # fm3 = fm4
 summary(fm4 <- update(fm2, random=list(ID=pdDiag(~TIME))))
 
 # 切片と傾き、相関あり
 summary(fm5 <- update(fm2, random=list(ID=pdSymm(~TIME))))
 
 # gls
 
 dose <- c(1,2,4,8,16)
 data.power <- data.frame(
   ID=rep(seq(6), each=5),
   dose=rep(dose, 6),
   cl=rep(1*exp(rnorm(6, sd=0.2)), each=5)
 )
 data.power$auc <- with(data.power,
   dose*1000/cl) * exp(rnorm(nrow(data.power), sd=0.2))
 
 data.power.g <- groupedData(auc ~ dose | ID, data.power)
 
 (fm11 <- gls(auc ~ factor(dose) -1, data.power))
 
 # dose 毎に異なる分散を推定
 (fm12 <- update(fm11, weights=varIdent(form=~ 1 | dose)))
 
 # ID 内で compound symmetry の相関
 (fm13 <- update(fm11, corr=corCompSymm(form=~1 | ID)))
 
 # ID 内でフルの相関
 (fm14 <- update(fm11, corr=corSymm(form=~1 | ID)))
 
 # nlme
 
 n <- 50
 m <- 5
 
 ccr <- rnorm(n, mean=80, sd=20)
 bw  <- rnorm(n, mean=60, sd=10)
 
 cl <- (ccr / 1000 + .01) * exp(rnorm(n, sd=0.2))
 vd <- (bw / 60 + 1) * exp(rnorm(n, sd=0.2))
 ke <- cl / vd
 
 data.ppk <- data.frame(
   ID=rep(seq(n), each=m),
   TIME=runif(n*m, min=0.5, max=24),
   CCR=rep(ccr, each=m),
   BW =rep(bw, each=m),
   CL =rep(cl, each=m),
   V  =rep(vd, each=m)
 )
 
 iv1 <- function(TIME, lCL, lV, Dose=1) {
   CL <- exp(lCL)
   V <- exp(lV)
   Dose / V * exp(-CL/V * TIME)
 }
 
 data.ppk$PRED <- with(data.ppk, iv1(TIME, log(CL), log(V)))
 data.ppk$DV <- data.ppk$PRED * exp(rnorm(nrow(data.ppk), sd=0.1))
 data.ppk.g <- groupedData(DV ~ TIME | ID, data.ppk, order=F,
   labels=list(x="Time", y="Concentration"),
   units=list(x="(hr)", y="(ng/mL)")
 )
 
 # lCL と lV にフルの相関
 (fm21 <- nlme(log(DV) ~ log(iv1(TIME, lCL, lV)), data.ppk.g,
   fixed=list(lCL~1, lV~1),
   random=(lCL+lV~1),
   start=c(log(0.1), log(1))
 ))
 
 # lCL と lV 独立
 (fm22 <- update(fm21, random=pdDiag(form=list(lCL~1, lV~1))))
 
 # lCL = TH1 + TH2 * CCR
 # start での初期値指定では切片が先に来る
 (fm23 <- update(fm22,
   fixed=list(lCL~CCR+1, lV~1),
   start=c(log(0.1), 0, log(1))
 ))
 
 # lV = TH3 + TH4 * BW
 (fm24 <- update(fm23,
   fixed=list(lCL~CCR+1, lV~BW+1),
   start=c(log(0.1), 0, log(1), 0)
 ))
 
 # ID 単位で CS 構造の個体内相関を考慮
 (fm25 <- update(fm22, corr=corCompSymm(form=~1|ID)))
 
 # 分散が予測値の 2 乗に比例と仮定
 # 2 (=1 * 2) に固定している。
 # デフォルトで予測値の a 乗に比例
 (fm26 <- update(fm22, weights=varPower(fixed=1)))