はじめに

機械学習を勉強するとまずはじめに遭遇するのが線形回帰や多重回帰といった線形モデルかと思います。線形と聞くと一直線の関数というイメージを持ちがちですが、この枠組みで線形といっているのはあくまでもパラメータに関して線形であるだけで、場合によってはグニャグニャ曲がった曲線も線形モデルに含まれます。例えば$y=ax^2+bx+c$は一見非線形ですが、パラメータに関しては線形であるため線形モデルに含まれます。このように線形モデルといってもかなり柔軟な表現が可能なわけですが、パラメータに関して非線形にすることでさらに表現力をアップさせたモデルがあり、それらは非線形モデルと呼ばれます。

この非線形モデルの特徴として①外挿に強いこと②パラメータ数を抑えることができること③パラメータに意味を見出すことができること、などがあり、色々な分野で使われています。また、線形モデルを改良したバージョンとして線形混合モデルがあるように、非線形モデルに対しても非線形混合モデルというモデルがあり、こちらもデータ数が少ない場合などよく用いられています。さらに、近年はゲノム情報を非線形モデルに適用しようという動きもあり、ゲノム情報を分散共分散行列に持つ非線形混合モデルというものも存在します。

ところがそれらをプログラミング言語Rで解こうとすると、様々なパッケージが存在するためどれを使えばよいか迷ってしまいます。今回は①非線形モデル②非線形混合モデル③ゲノム情報を含む非線形混合モデルそれぞれに対し、Rで使えるパッケージを紹介＆比較をしていこうと思います。なお、③のゲノム情報を含む非線形混合モデルはかなり特殊な内容です。その分野でない方は②の非線形混合モデルまで読めば十分だと思います。なお、非線形混合モデルをベイズ的に解いたモデルがあります。厳密にはこれらは非線形混合モデルとは別ですが、ベイズ的に解いているパッケージも②、③の中に入れさせていただきました。

非線形モデル

まずは最も基本的な変量効果を含まない非線形モデルを解く際に使えるパッケージを比較していきます。

データ生成

以下ではいくつかのパッケージを用いながらデータへのフィッティングを行いますが、その際用いるデータは簡単のため、3パラメータロジスティック曲線 $y=\frac{A}{1+B\exp(-Ct)}$を用いて人工的に生成したものを使いました。生成は以下のように行っています。

# それぞれのパラメータの値を指定
A <- 10.0
B <- 10.0
C <- 0.5

# データ生成
time <- seq( 0.0, 10.0, 0.1 )
y_true <- A / ( 1 + B*exp( -C*time ) )　# ノイズを含まない真値
y_real <- y_true + rnorm( length( time ), 0.0, 1.0 ) #　ノイズを含む値

# プロット
plot( time, y_true )
plot( time, y_real )

⇓こちらが真値のプロット

⇓こちらがノイズを含んだデータのプロット

ちょっとノイズかけすぎている感じも否めないのですがいったんはこのデータを使ってフィッティングを行っていきます。
では実際にパラメータを比較していきましょう。

nls

非線形を扱うパッケージの中でおそらく最も有名なのがこのnls関数。ベース関数の中に入っているので特にパッケージをインストールすることなく使うことができます。
引数として必ず必要なのは以下の3つです。
①式
②データ(基本的にデータフレーム)
③初期値

実際に使いながらその使い心地を見ていきましょう。

# データフレームの用意
data_nl <- data.frame( x = time, y = y_real )

# 結果の格納 (初期値はそれぞれA=1,B=1,C=1で指定)
result_nl <- nls( formula = y ~ a / ( 1 + b*exp( -c*x ) ), 
                  start = c( a = 1, b = 1, c = 1),
                  data = data_nl )

# 結果の表示
summary( result_nl )

推定された結果はこんな感じ

Formula: y ~ a/(1 + b * exp(-c * x))

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a  9.58438    0.37742  25.394  < 2e-16 ***
b 10.16359    1.68520   6.031 2.89e-08 ***
c  0.55432    0.05294  10.472  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9973 on 98 degrees of freedom

Number of iterations to convergence: 8 
Achieved convergence tolerance: 1.301e-06

なかなか良い精度で推定ができています。なお、nlsはデータフレームを指定せず次のように直接ベクトルを指定することも可能です。

result_nl <- nls( formula = y_real ~ a / ( 1 + b*exp( -c*time ) ), 
                  start = c( a = 1, b = 1, c = 1) )

さて、このnlsでは初期値をA=1,B=1,C=1と指定しましたが、うまく良い値に指定できないと推定ができずエラーを吐いてしまいます。今回は人工的に生成したデータということもあり結構離れた値でもうまく推定ができましたが生データだともう少しシビアになる印象です。私はGeoGebraでフィッティングしながらよさげな値を探していましたが、結構骨の折れる作業です。そんな初期値の設定という問題からの解法を目指して作られたのが次のdrm関数です。

drm

drcパッケージというパッケージの中に入っている関数です。nlsと違い初期値を必要としないという特徴があります。なお、初期値を必要としないとは言ったものの決して初期値を使っていないわけではありません。データの特徴から初期値を推定して使っています。詳しくは次のサイトをご覧ください。
https://www.r-bloggers.com/2020/02/a-collection-of-self-starters-for-nonlinear-regression-in-r/

この関数のもう一つのありがたいところはnlsのように非線形の式をわざわざ手で与える必要がないという点です。あらかじめ使える非線形モデルに名前が付けられており、fctという引数にその名前を指定すればその関数にフィッティングしてくれます。3パラメータのロジスティック曲線はL.3()という名前がついており、次のように書けば推定してくれます。

# 結果の格納
result_drm <- drm( data = data_nl,
                   fct = L.3() )

# 結果の表示
summary(result_drm)

推定された結果はこんな感じ

Model fitted: Logistic (ED50 as parameter) with lower limit fixed at 0 (3 parms)

Parameter estimates:

               Estimate Std. Error t-value   p-value    
b:(Intercept) -0.405881   0.043771 -9.2729 4.512e-15 ***
d:(Intercept) 11.063374   0.964862 11.4663 < 2.2e-16 ***
e:(Intercept)  6.184889   0.529027 11.6911 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error:

 0.9930205 (98 degrees of freedom)

ちょっとわかりにくいのですが、-bが先ほどのC、dが先ほどのA、eが先ほどのBに対応しています。nlsよりは予測精度が下がってしまいましたが、初期値を指定しなくても推定ができることが分かるかと思います。

非線形混合モデル

続いて非線形混合モデルです。いくつかの個体に対して時系列データを持っている場合、それぞれに非線形モデルでフィッティングを行うと、個体ごとに異なるパラメータの値$A,B,C$が得られるかと思います。個体iのパラメータを$A_{i},B_{i},C_{i}$としたとき、個体数が増えた場合、この$A_{i},B_{i},C_{i}$は何らかの分布に従うことが考えられます。そしてその分布として正規分布を仮定したのが非線形混合モデルと呼ばれるモデルです。つまり$A_{i}\overset{}\sim N(\mu_A,\sigma_A^2),B_{i}\sim N(\mu_B,\sigma_B^2),C_{i} \sim N(\mu_C,\sigma_C^2)$を仮定しています。

数学的な定式化は今回の趣旨からずれるのでこの程度にしますが、このように全個体が正規分布により結びつくため、データへのフィッティングは全個体同時に行います。従ってそれに合うようにデータフレームを作る必要があることに注意が必要です。

なお、以下で紹介する内容は次の論文で紹介されている物をベースにしています。必要があればこちらをご覧ください。
https://www.tandfonline.com/doi/full/10.1080/10705511.2017.1396187

データ生成

まずデータを生成します。非線形混合モデルではパラメータを確率変数と見なして推定するため一個体ではなかなか推定が困難です。そこで今回は10個体分のデータを人工的に生成しました。

library( mvnfast )
library( tidyverse )

# mean for A, B, C
A <- 10.0
B <- 10.0
C <- 1.0

# variance for A, B, C
var_A <- 1.0
var_B <- 1.0
var_C <- 0.1

# create vectors for means and variances
mu <- c( A, B, C )
S<-diag( c( var_A, var_B, var_C ) )

# number of individuals
N <- 10

# generate parameters for each individual
Par <- mvnfast::rmvn( N, mu, S )

# define logistic model to create phenotype
Logistic <- function( i, x, par ){
  return( par[ i, 1 ] / ( 1 + par[ i, 2 ] * exp( - par[ i, 3 ] * x ) ) )
}

# time sequence
time <- seq( 0.0, 10.0, 0.1 )

# generate phenotype without random error
u <- matrix( unlist( map( 1 : N, Logistic, time, Par) ), ncol = length( time ), byrow = TRUE )

# generate random error
e <- matrix( rnorm( N * length( time ), 0, 0.1 ), nrow = N )

# phenotype with random error 
y <- u + e

matplot( time, t( u ), type = "l", ylab = "phenotype" )
matplot( time, t( y ), type = "l", ylab = "phenotype" )

⇓こちらが生成した真の曲線

⇓こちらがノイズを乗せた曲線

データはyという10(個体数)×101(時系列観測数)の行列に格納されています。ではこのデータを実際に用いて推定を行っていきます。

nlme

まずはじめにnlmeというパッケージのnlme関数を用いてフィッティングを行います。
nlsと同様nlmeでも①式②データ(データフレーム)③初期値が必要になりますが非線形混合モデルではそれに加えてfixed,randomという引数が出てきます。
非線形混合モデルでは、例えばパラメータ$A_i$は$A_i=\mu_A+a_i$というように平均$\mu_A$とそこからのずれ$a_i$の和でかかれます。$a_i$は平均0,分散$\sigma_g^2$に従う確率変数です。この時$\mu_A$は固定効果(fixed effect)、$a_i$は変量効果(random effect)と呼ばれます。今回はパラメータが$A_i,B_i,C_i$の3つありますが、例えば$B_i$については変量効果だけ考えて固定効果を考えない($B_i=0+b_i$)、なんてこともできるわけです。そういうのを指定するのがfixed,randomという引数です。今回はA,B,Cの3つ全てがfix,randomに指定されていますので、すべてのパラメータで固定効果、変量効果両方が考慮されていることになります。
なお、y_nlmeは以下のような形のデータフレームです。

> head(y_nlme)
         y   t    id
1 1.068608 0.0 ind 1
2 1.177144 0.1 ind 1
3 1.285698 0.2 ind 1
4 1.319774 0.3 ind 1
5 1.570552 0.4 ind 1
6 1.486793 0.5 ind 1

以下はnlmeを用いたコード例です。データフレームをyから作り出し、nlmeを用いてフィッティングを行っています。また初期値は本来nls等を用いて推定するべきですが、今回は真値に近い値を使わせてもらいました。結構、初期値依存性が強く、ちょっとでも真値から離れると推定できなくなってしまうようです。

ibrary(nlme)

y_nlme <- data.frame( y = as.vector( t( y ) ),
                      t = rep( time, N ),
                      id = rep( paste( "ind", 1:N ), each = length( time ) )
                      )

result_nlme <- nlme( model = y ~ A / ( 1 + B * exp( - C * t ) ),
                     data = y_nlme,
                     fixed = A + B + C ~ 1,
                     random = A + B + C ~ 1,
                     group = ~ id,
                     start = c( 9, 9, 1.2) )

summary( result_nlme )

結果はこんな感じです↓

Nonlinear mixed-effects model fit by maximum likelihood
  Model: y ~ A/(1 + B * exp(-C * t)) 
  Data: y_nlme 
        AIC       BIC   logLik
  -1632.555 -1583.378 826.2776

Random effects:
 Formula: list(A ~ 1, B ~ 1, C ~ 1)
 Level: id
 Structure: General positive-definite, Log-Cholesky parametrization
         StdDev     Corr         
A        1.06282666 A      B     
B        0.72486442  0.372       
C        0.38344919  0.021 -0.636
Residual 0.09671878              

Fixed effects:  A + B + C ~ 1 
      Value Std.Error  DF  t-value p-value
A 10.012274 0.3369050 998 29.71839       0
B  9.913340 0.2368246 998 41.85942       0
C  0.953666 0.1214720 998  7.85091       0
 Correlation: 
  A      B     
B  0.359       
C  0.020 -0.611

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-3.20176694 -0.65692881 -0.03407789  0.67654217  3.27744939 

Number of Observations: 1010
Number of Groups: 10

$(\mu_A,\mu_B,\mu_C)=(10,10,1)$、$(\sigma_A,\sigma_B,\sigma_C)=(1,1,0.1)$と設定していて推定値が$(\mu_A,\mu_B,\mu_C)=(10.0,9.9,0.95)$、$(\sigma_A,\sigma_B,\sigma_C)=(1.1,0.7,0.4)$なのでなかなか良い推定と言えるのではないでしょうか

なお、推定された各個体のパラメータ$A_i,B_i,C_i$はresult_nlme$coefficientsに格納されています。中身を見てみるとこんな感じ↓

> result_nlme$coefficients
$fixed
         A          B          C 
10.0122743  9.9133396  0.9536662 

$random
$random$id
                A           B          C
ind 1   0.4170000 -0.20138288  0.1669172
ind 10  0.6099252  0.89173653 -0.3487278
ind 2  -2.0442634 -0.49581776 -0.3926657
ind 3   1.1958877  0.09492872 -0.4991427
ind 4  -1.5691441 -0.98687613  0.3428136
ind 5   0.6228087 -0.54745113  0.4318851
ind 6   1.1732598 -0.43455494  0.4601597
ind 7   0.2829085  1.52358578 -0.5623956
ind 8  -0.8811667 -0.11911876  0.1545433
ind 9   0.1927842  0.27495058  0.2466129

推定された$A_i$と真の$A_i$を比べてみるとこんな感じ。かなり良く推定できていることが分かるかと思います。

nlmer(lme4)

次にlme4パラメータに入っているnlmer関数を用いて推定を行います。nlmeがテイラー展開による一次の項までの近似でアルゴリズムを回しているのに対しnlmerではラプラス近似と呼ばれる方法によりアルゴリズムを回しているという違いがあります。
nlmeと同様①式②データ(データフレーム)③初期値の3つが必要になりますが、nlmeと違い式に少し工夫をする必要があります。工夫というのはderiv関数という関数にかけて微分した式も準備しておく点です。以下のコードではlog_modelが元の式、log_modelgがderiv関数にかけて微分もできるようにした式です。~ ( A + B + C | id )の部分は変量効果を指定しており|idはidごとにパラメータA,B,Cの値が違うことを表しています。

library( lme4 )

log_model <- function( t, A, B, C )  A / ( 1 + B * exp( - C * t ) )

log_modelg <- deriv( body( log_model ),
                       namevec = c( "A", "B", "C" ),
                       function.arg = log_model )

result_nlmer <- nlmer( y ~ log_modelg( t, A, B, C ) ~ ( A + B + C | id ),
                       data = y_nlme,
                       start = c( A = 9, B = 9, C = 1.5 ) )

summary( result_nlmer )

⇓こちらが計算結果のまとめです。

Nonlinear mixed model fit by maximum likelihood  ['nlmerMod']
Formula: y ~ log_modelg(t, A, B, C) ~ (A + B + C | id)
   Data: y_nlme

     AIC      BIC   logLik deviance df.resid 
 -1378.2  -1329.0    699.1  -1398.2     1000 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.6403 -0.6675  0.0262  0.6650  3.6410 

Random effects:
 Groups   Name Variance Std.Dev. Corr     
 id       A    0.88838  0.9425            
          B    0.02978  0.1726   1.00     
          C    0.07122  0.2669   0.43 0.43
 Residual      0.01246  0.1116            
Number of obs: 1010, groups:  id, 10

Fixed effects:
  Estimate Std. Error t value
A  9.77622    0.29823   32.78
B 10.03667    0.08667  115.80
C  0.91133    0.08444   10.79

Correlation of Fixed Effects:
  A     B    
B 0.622      
C 0.426 0.291
警告メッセージ:
1: vcov.merMod(object, use.hessian = use.hessian) で:
  variance-covariance matrix computed from finite-difference Hessian is
not positive definite or contains NA values: falling back to var-cov estimated from RX
2: vcov.merMod(object, correlation = correlation, sigm = sig) で:
  variance-covariance matrix computed from finite-difference Hessian is
not positive definite or contains NA values: falling back to var-cov estimated from RX

$B_i$の標準偏差$\sigma_B$の推定精度がnlmeより悪くなってしまった気がします。警告も出ていますね。この警告はヘッシアンに基づいて計算した分散共分散行列が半整定値にならなかったためRXと呼ばれれる別の行列を代用したことを言っています。これが推定精度の低下に影響を与えたかどうかは分かりませんが、ちょっと気になりますね。

saemix

次にsaemixパッケージの中のsaemix関数です。この関数はEMアルゴリズムを改良したstochastic approximattion EMアルゴリズム(SAEM)を用いています。例えば式やデータをsaemixパッケージが適用するオブジェクトに変換しないといけないなど、nlmeなどと比べて準備や引数が多くややこしいのですが、その分様々な機能が提供されており使い慣れればかなり便利だろうと思います。
今までは①式②データ(データフレーム)③初期値を準備していましたが、今回は①saemix formatに沿って準備されたモデル②SaemixData object形式のデータの2つを準備する必要があります。なお、初期値が不要なわけではなく①の中に含まれています。
以下のコードではjbが②に対応し、saemix.modelが①に対応しています。jbはsaemixDataという関数を用いて用意します。その関数の引数はデータフレームの名前、グループを指定している列名、独立変数の列名、従属変数の列名です。saemix.modelでは式、初期値の行列、各パラメータの分布(transform.par)、パラメータ推定の有無(fixed.estim)等を引数として渡します。式は別途定義されている必要があり、それはjb.modelで指定しています。ちょっと変わった式定義ですが、この引数などは変えてはいけないものと思われます。
また、最後にoptというオブジェクトを定義していますが、これはランダム性に関してシードを固定したり、MCMCの記録を残すかどうかなどを決めているオプションです。これらすべてを合わせてsaemixという関数に突っ込めば推定が始まります。

library(saemix)

jb <-saemixData( name.data = y_nlme,
                 name.group = "id",
                 name.predictors = "t",
                 name.response = "y" )

jb.model <-function( psi, id, xidep ){
  time <- xidep[ , 1 ]
  A <- psi[ id, 1 ]
  B <- psi[ id, 2 ]
  C <- psi[ id, 3 ]
  resp <- A / ( 1 + B * exp( -C * time ) )
  return( resp )
}

saemix.model <- saemixModel( model = jb.model,
                             psi0 = matrix( c( 9, 9, 1.5 ), ncol = 3, byrow = TRUE,
                                            dimnames = list( NULL, c( "A", "B", "C" ) ) ),
                             transform.par = c( 0, 0, 0 ),
                             fixed.estim = c( 1, 1, 1 ) )

opt <- list(seed = 1, save = FALSE, save.graphs = FALSE )

result.jb <- saemix( saemix.model, jb, opt )

print( result.jb )

こちらが結果です↓

Nonlinear mixed-effects model fit by the SAEM algorithm
-----------------------------------
----          Data             ----
-----------------------------------
Object of class SaemixData
    longitudinal data for use with the SAEM algorithm
Dataset y_nlme 
    Structured data: y ~ t | id 
    Predictor: t () 
Dataset characteristics:
    number of subjects:     10 
    number of observations: 1010 
    average/min/max nb obs: 101.00  /  101  /  101 
First 10 lines of data:
      id   t         y mdv cens occ ytype
1  ind 1 0.0 0.8287020   0    0   1     1
2  ind 1 0.1 0.7439885   0    0   1     1
3  ind 1 0.2 0.8188963   0    0   1     1
4  ind 1 0.3 0.9632055   0    0   1     1
5  ind 1 0.4 0.9861669   0    0   1     1
6  ind 1 0.5 1.1284166   0    0   1     1
7  ind 1 0.6 1.1208325   0    0   1     1
8  ind 1 0.7 1.3280937   0    0   1     1
9  ind 1 0.8 1.3193079   0    0   1     1
10 ind 1 0.9 1.5355938   0    0   1     1
-----------------------------------
----          Model            ----
-----------------------------------
Nonlinear mixed-effects model
  Model function
  Model type:  structural
function( psi, id, xidep ){
  time <- xidep[ , 1 ]
  A <- psi[ id, 1 ]
  B <- psi[ id, 2 ]
  C <- psi[ id, 3 ]
  resp <- A / ( 1 + B * exp( -C * time ) )
  return( resp )
}
<bytecode: 0x5587026f8630>
  Nb of parameters: 3 
      parameter names:  A B C 
      distribution:
     Parameter Distribution Estimated
[1,] A         normal       Estimated
[2,] B         normal       Estimated
[3,] C         normal       Estimated
  Variance-covariance matrix:
  A B C
A 1 0 0
B 0 1 0
C 0 0 1
  Error model: constant , initial values: a.1=1 
    No covariate in the model.
    Initial values
             A B   C
Pop.CondInit 9 9 1.5
-----------------------------------
----    Key algorithm options  ----
-----------------------------------
    Estimation of individual parameters (MAP)
    Estimation of standard errors and linearised log-likelihood
    Estimation of log-likelihood by importance sampling
    Number of iterations:  K1=300, K2=100 
    Number of chains:  5 
    Seed:  1 
    Number of MCMC iterations for IS:  5000 
    Simulations:
        nb of simulated datasets used for npde:  1000 
        nb of simulated datasets used for VPC:  100 
    Input/output
        save the results to a file:  FALSE 
        save the graphs to files:  FALSE 
----------------------------------------------------
----                  Results                   ----
----------------------------------------------------
-----------------  Fixed effects  ------------------
----------------------------------------------------
     Parameter Estimate SE     CV(%)
[1,] A          9.782   0.2994 3.1  
[2,] B         10.128   0.3231 3.2  
[3,] C          0.913   0.0853 9.3  
[4,] a.1        0.099   0.0022 2.2  
----------------------------------------------------
-----------  Variance of random effects  -----------
----------------------------------------------------
  Parameter Estimate SE    CV(%)
A omega2.A  0.896    0.401 45   
B omega2.B  1.005    0.467 46   
C omega2.C  0.073    0.033 45   
----------------------------------------------------
------  Correlation matrix of random effects  ------
----------------------------------------------------
         omega2.A omega2.B omega2.C
omega2.A 1        0        0       
omega2.B 0        1        0       
omega2.C 0        0        1       
----------------------------------------------------
---------------  Statistical criteria  -------------
----------------------------------------------------
Likelihood computed by linearisation
      -2LL= -1619.158 
      AIC = -1605.158 
      BIC = -1603.04 

Likelihood computed by importance sampling
      -2LL= -1605.372 
      AIC = -1591.372 
      BIC = -1589.254 
----------------------------------------------------

色々書いてありますがResultsのところを見ると推定値が書いてありました。

----------------------------------------------------
----                  Results                   ----
----------------------------------------------------
-----------------  Fixed effects  ------------------
----------------------------------------------------
     Parameter Estimate SE     CV(%)
[1,] A          9.782   0.2994 3.1  
[2,] B         10.128   0.3231 3.2  
[3,] C          0.913   0.0853 9.3  
[4,] a.1        0.099   0.0022 2.2  
----------------------------------------------------
-----------  Variance of random effects  -----------
----------------------------------------------------
  Parameter Estimate SE    CV(%)
A omega2.A  0.896    0.401 45   
B omega2.B  1.005    0.467 46   
C omega2.C  0.073    0.033 45

確かにいい感じに推定されていますね。パラメータ間の共分散なども定義できるようなのでnlmeやnlmerよりはもう少し柔軟なモデル設計ができそうです。

brms

最後に紹介するのがベイジアン非線形モデルを持ちいたパッケージbrmsのbrms関数です。ベイズ統計によるモデリングはあまり詳しくないのですが特に共役事前分布などは仮定せず自由に事前分布を設定しメトロポリスヘイスティングによるMCMCを回しているんだと思います。
コードは比較的シンプルで、①変量効果を指定した式②データフレーム③事前分布の定義の3つさえあれば動きます。
おそらく見ればすぐ分けると思いますが、唯一言及すべきは変量効果の指定方法でしょう。list( A ~ (1|2|id), B ~ (1|2|id), C ~ (1|2|id) )というコードがありますがここが変量効果の指定をしています。(1|2|id)の部分が対象のパラメータが変量効果を含みidごとに違う値を取ることを示しています。

library( brms )

log_brms <- y ~ A / ( 1 + B * exp( - C * t ) )

form = bf( log_brms, nl=TRUE ) + list( A ~ (1|2|id), B ~ (1|2|id), C ~ (1|2|id) )

prior <- c( set_prior( "normal(9, 1.5)", nlpar = "A" ),
            set_prior( "normal(9, 1.5)", nlpar = "B" ),
            set_prior( "normal(1.5, 1)", nlpar = "C" ) )

result_brms <- brm( form, data = y_nlme, prior = prior )

summary( result_brms )

回している間には次のようなメッセージが出力されます。

Compiling Stan program...
Start sampling

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000469 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 4.69 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 271.942 seconds (Warm-up)
Chain 1:                316.819 seconds (Sampling)
Chain 1:                588.761 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000317 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 3.17 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 119.35 seconds (Warm-up)
Chain 2:                198.953 seconds (Sampling)
Chain 2:                318.303 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000318 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 3.18 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 89.719 seconds (Warm-up)
Chain 3:                226.889 seconds (Sampling)
Chain 3:                316.608 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000398 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 3.98 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 105.094 seconds (Warm-up)
Chain 4:                318.744 seconds (Sampling)
Chain 4:                423.838 seconds (Total)
Chain 4: 
警告メッセージ:
1: There were 834 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them. 
2: There were 3017 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded 
3: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See
https://mc-stan.org/misc/warnings.html#bfmi-low 
4: Examine the pairs() plot to diagnose sampling problems
 
5: The largest R-hat is 4.95, indicating chains have not mixed.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#r-hat 
6: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess 
7: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

ただこの関数、ベイズ的に実装しているだけあって推定にスーパー時間がかかります。事前分布を自由に設定できるなど、非線形混合モデルにはない点はあるので、うまくすみ分けるのが大切なように思います。さて、下はその推定結果です。

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: y ~ A/(1 + B * exp(-C * t)) 
         A ~ (1 | 2 | id)
         B ~ (1 | 2 | id)
         C ~ (1 | 2 | id)
   Data: y_nlme (Number of observations: 1010) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Multilevel Hyperparameters:
~id (Number of levels: 10) 
                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(A_Intercept)                322.02    829.50     0.01  2649.17 4.77        4        4
sd(B_Intercept)               1817.96   5741.47     0.22 23062.55 4.96        4       11
sd(C_Intercept)                  1.33      1.02     0.25     2.96 3.49        4       NA
cor(A_Intercept,B_Intercept)    -0.45      0.53    -0.95     0.54 2.53        5        4
cor(A_Intercept,C_Intercept)     0.25      0.25    -0.34     0.63 1.98        9       35
cor(B_Intercept,C_Intercept)    -0.22      0.54    -0.98     0.68 2.68        5       11

Regression Coefficients:
            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
A_Intercept     3.24      3.93    -0.24    10.20 3.53        5       24
B_Intercept     4.07      5.39    -2.31    10.90 2.47        5        4
C_Intercept     0.18      1.46    -2.16     1.91 4.03        4       NA

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     4.72      2.98     0.10     7.91 3.36        5       29

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
警告メッセージ:
1: Parts of the model have not converged (some Rhats are > 1.05). Be careful when analysing the results! We recommend running more iterations and/or setting stronger priors. 
2: There were 834 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

どうも推定値が変だと思って色々考えていましたが、そもそも事前分布の設定等色々違う気がしてきました…どうもベイジアンへの理解がまだ微妙で…いずれ正しいものに更新するので一旦はこんな感じで失礼しますm(__)m

非線形混合 (ゲノム関係行列の指定アリ)

さて、個体ごとに得られる3つのパラメータ$A_i,B_i,C_i$について、遺伝的に近い個体同士は近しい値を取ると考えられます。これは遺伝的な近さを表す行列を$K$としたとき$\mathbf{A}=(A_1,A_2,..A_n),\mathbf{B}=(B_1,B_2,..B_n),\mathbf{C}=(C_1,C_2,..C_n)^T$がそれぞれ$\mathbf{A}\sim MVN(\mathbf{1}\mu_A,K\sigma_A^2),\mathbf{B}\sim MVN(\mathbf{1}\mu_B,K\sigma_B^2),\mathbf{C}\sim MVN(\mathbf{1}\mu_C,K\sigma_C^2)$に従うことと同値です。この遺伝的な近さを表す行列$K$はゲノム関係行列と呼ばれゲノム情報から計算可能です。NGS(次世代シーケンサー)の発達によりゲノム情報がより簡単に手に入るようになった今、このようなモデルは育種(品種改良)において重要になってきています。ここではそのようなモデル特化したパッケージとしてGenomeBasedModelとSAEMによるパッケージの二つを紹介します。

データ生成

まずデータを生成しますがこれはほぼ非線形混合モデルを紹介した時と同じです。ゲノム関係行列はここで計算すると大変なので、約1500個体のイネの関係行列をGRMというオブジェクトに入れておきました。突然出てきますがご了承ください。

library( mvnfast )
library( tidyverse )

# mean for A, B, C
A <- 10.0
B <- 10.0
C <- 1.0

# variance for A, B, C
var_A <- 1.0
var_B <- 1.0
var_C <- 0.1

# number of individuals
N <- 10

# genomic relationship matrix of N individuals
K <- GRM[ 1:N, 1:N ]

# generate parameters for each individual
u_a <- mvnfast::rmvn( 1, rep( A, N ), K * var_A )
u_b <- mvnfast::rmvn( 1, rep( B, N ), K * var_B )
u_c <- mvnfast::rmvn( 1, rep( C, N ), K * var_C )

# define logistic model to create phenotype
Logistic <- function( i, x, par ){
  return( par[ i, 1 ] / ( 1 + par[ i, 2 ] * exp( - par[ i, 3 ] * x ) ) )
}

# time sequence
time <- seq( 0.0, 10.0, 0.1 )

# generate phenotype without random error
u <- matrix( unlist( map( 1 : N, Logistic, time, Par) ), ncol = length( time ), byrow = TRUE )

# generate random error
e <- matrix( rnorm( N * length( time ), 0, 0.1 ), nrow = N )

# phenotype with random error 
y <- u + e

matplot( time, t( u ), type = "l", ylab = "phenotype" )
matplot( time, t( y ), type = "l", ylab = "phenotype" )

以下の写真はそれぞれノイズを含まない真の曲線とノイズを含む表現型データです。挙動は非線形混合モデルの時とほぼ同じですね。

GenomeBasedModel

GenomeBasedModelはこちらのサイトで公開されているRパッケージGenomeBasedModelの中にある関数です。変分ベイズとMCMCを組み合わせることにより推定を行っています。

library( GenomeBasedModel )

################################################################################
# prepare the dataset for estimation

Input <- matrix( time, nrow = length( time ), ncol = N ) # independent variable for each individual
Y <- rbind( 1 : N, t( y ) ) # dependent variable for each individual
Geno <- rbind( 1 : N, rep( 1, N ) ) # not important but necessary
Kind <- rbind( 1 : N, K ) # genomic relationship matrix
missing.value <- 9999 # number indicating missin value
Init <- c( 9, 9, 1.5 ) # initial value

################################################################################
# definition of logistic function for GenomeBased Model

GBM_logit <- function( input, freevec, parameters ) {
  A <- parameters[ 1 ]
  B <- parameters[ 2 ]
  C <- parameters[ 3 ]
  y <- A / ( 1 + B * exp( - C * input ) )
  return( y )
}

################################################################################
# parameter estimation

result <- GenomeBasedModel( Input = Input, Freevec = 0, Y = Y,
                            Missing = missing.value, Np = 3, Geno = Geno,
                            Methodcode = c( 8, 8, 8 ), Referencevalues = Init,
                            Modelfunction = GBM_logit, K=Kind )
# calculate means of parameters for each individual from posterior distribution
A.est <- apply( result$Para[[ 1 ]], 2, mean)
B.est <- apply( result$Para[[ 2 ]], 2, mean)
C.est <- apply( result$Para[[ 3 ]], 2, mean)

# calculate global means of each parameter
A.mean <- mean( A.est )
B.mean <- mean( B.est )
C.mean <- mean( C.est )

さて、推定されたパラメータを出力した結果がこちらです↓

> print( A.mean )
[1] 9.900656
> print( B.mean )
[1] 10.65743
> print( C.mean )
[1] 1.165277

いい感じに推定されています。推定された値と真値をプロットしてみてもびしっと$y=x$上に乗っていることが分かりますね

SAEM (Maud)

こちらも紹介できればと思ったのですが、パッケージ化されていないこともあり、コードが長くなりすぎるので断念しました。詳しくはこちらの論文で紹介されているので興味がある方は読んでみてください。関数はこちらから手に入れることができます。