TMBで、打ちきりデータへのモデルのあてはめ

TMBを使って、打ちきりデータへのモデルのあてはめをやってみます。Stan User’s GuideのCensored Dataの項の内容をTMBに移植したような内容になっています。

データ

模擬データを生成します。平均10、標準偏差4の正規分布で、6以下の値が打ち切られて欠損値になっているという状況を想定しています。

set.seed(123)
N <- 240
mu <- 10
sigma <- 4
cens <- 6
Y <- rnorm(N, mu, sigma)

平均と標準偏差を見てみます。

mean(Y)
## [1] 9.978479
sd(Y)
## [1] 3.753209

打ちきり閾値以下の値を欠損値にしたものをY2とします。

Y2 <- ifelse(Y <= cens, NA, Y)

欠損値の数は以下のようになっています。

length(Y[is.na(Y2)])
## [1] 37

打ちきりを考慮せず、単純に欠損値を除いた平均と標準偏差の値です。打ちきりにより、平均は大きく、標準偏差は小さくなっています。

mean(Y2, na.rm = TRUE)
## [1] 10.95027
sd(Y2, na.rm = TRUE)
## [1] 3.198584

ヒストグラム

ヒストグラムを見てみます。青緑が打ち切られた部分になります。

library(ggplot2)

data.frame(Y = Y, censored = (Y <= cens)) |>
  ggplot(aes(x = Y, fill = censored)) +
  geom_histogram(binwidth = 1, center = 0.5, na.rm = TRUE) +
  theme_minimal()

モデル

TMBのC++のコードです。"models/censored.cpp"に保存しておきます。TMBのpnorm関数は、Rのpnormと同様に、正規分布の累積分布関数です。TMBのdnormには、対数で値を返すかどうかを指定するgive_log引数があるのですが、pnormにはありません（ドキュメント）。試してみたところ、対数で値を返すようになっていました。

censored.cpp

// Censored data
#include <TMB.hpp>

template<class Type>
Type objective_function<Type>::operator() ()
{
  DATA_VECTOR(Y_obs);
  DATA_SCALAR(L);
  DATA_SCALAR(N_cens);
  PARAMETER(mu);
  PARAMETER(sigma);
  Type nll = 0;

  nll += -sum(dnorm(Y_obs, mu, sigma, true));
  nll += -N_cens * pnorm(L, mu, sigma);
  return nll;
}

コンパイルと最適化

モデルをコンパイルして、できたライブラリをロードします。

library(TMB)

model_name <- "censored"
file.path("models", paste(model_name, "cpp", sep = ".")) |>
  compile()
file.path("models", dynlib(model_name)) |>
  dyn.load()

MakeADFun関数で、最適化関数に渡すオブジェクトを作成して、nlminb関数で最適化します。

data <- list(Y_obs = Y2[!is.na(Y2)], L = cens,
             N_cens = N - length(Y2[is.na(Y2)]))
parameters <- list(mu = 1, sigma = 1)
obj <- MakeADFun(data, parameters, DLL = model_name)
## Constructing atomic pnorm1
opt <- nlminb(obj$par, obj$fn, obj$gr)
## outer mgc:  21962.24 
## outer mgc:  2635.305 
## outer mgc:  1383.968 
## outer mgc:  650.5126 
## outer mgc:  444.6701 
## outer mgc:  157.9629 
## outer mgc:  90.96684 
## outer mgc:  44.74458 
## outer mgc:  22.20723 
## outer mgc:  9.387689 
## outer mgc:  3.270364 
## outer mgc:  2.659021 
## outer mgc:  1.760704 
## outer mgc:  1.038897 
## outer mgc:  0.4863183 
## outer mgc:  0.101038 
## outer mgc:  0.01264069 
## outer mgc:  0.001014593 
## outer mgc:  4.565342e-05

結果

結果です。convergenceが0なので、収束しているようです。muとsigmaは、打ちきりを考慮しない場合よりも、元の値に近くなっていました。

print(opt)
## $par
##        mu     sigma 
## 10.114447  3.783921 
## 
## $objective
## [1] 505.7071
## 
## $convergence
## [1] 0
## 
## $iterations
## [1] 18
## 
## $evaluations
## function gradient 
##       23       19 
## 
## $message
## [1] "relative convergence (4)"