lmeとnlme - Archives@Wiki

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

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)))