procでArrayステートメントを使わずにDoループを使う方法SASのNLMIXED

Question

次のsasコードを実行すると、最適化が完了できません。問題の1つは、私は推測することができますが、私はproc NLMIXEDには適していないワイドフォーマットのデータを使用している、誰も私は、次のコード/問題、しかし、どのように私は長い形式のデータでそれを行うことができます。

proc nlmixed 
    data=new method=GAUSS NOAD qpoints=50 tech=newrap maxit=5000 
; 
    array y{*} y1-y5; 
    array Time{*} Time1-Time5; 
    parms b1=0.95 sdn1=0.95 sdC=0.90 fi=0.45; 
    bounds sdn1 >= 0, sdC >= 0; 
    pi=constant('pi'); 
    L1=log(pi)-((pi*(y[1]-b1*Time[1] + bD))/(sdn1*sqrt(3)))-log(sdn1)-0.5*log(3)-2*log(1+exp(-((pi*(y[1]-b1*Time[1] + bD))/(sdn1*sqrt(3))))); 
    L2=1; 
    do j=2 to 5; 
    L2=L2*(log(pi)-((pi*(y[j]-b1*Time[j]+y[j-1]*fi+ bD))/(sdn1*sqrt(3*(1-fi*fi))))-log(sdn1)-0.5*log(3*(1-fi*fi))-2*log(1+exp(-((pi*(y[j]-b1*Time[j]+y[j-1]*fi+ bD))/(sdn1*sqrt(3*(1-fi*fi))))))); 
    end; 
    L=L1*L2;dum=1; 
    model dum ~ general(L); 
    random bD ~ normal(0,sdC*sdC) subject=ID; 
run; 

Answer 1

私はポイントのいくつかを持っています。まず、複雑な尤度関数を使用して問題が発生していると思われる場合は、データステップで尤度関数を作成することによって問題の原因を突き止めることができます。オーバーフローの原因となる値を指数化しようとしているか、または負の数の対数または平方根を計算しようとする状況があります。エラー報告は、データステップコードで（ログ）尤度を構築しようとすると、そのような問題についてより多くの情報を提供します。データステップコードを書く目的は、可能性を最大限にすることではなく、指定されたパラメータセットの可能性を評価することです。

NLMIXEDコードをデータステップコードに変換するためにいくつかの変更が必要です。まず、NLMIXEDステートメント（オプション付き）をDATAステートメントに変換する必要があります。 AデータNULL;通常は十分です。次に、PARMSステートメントはDATAステップではサポートされていません。これをRETAIN文に変換することができます。これはしばしばコードを変更します。 PARMSステートメントでは、複数の値の指定をサポートしません（PARMSステートメントでは、パラメーターと値の間で等号が許されます）他の文は、DATAステップでは使用できないNLMIXEDコードでよく発生します。あなたのコードでは、BOUNDS、MODEL、およびRANDOMステートメントをコメントアウトする必要があります。残念なことに、RANDOMステートメントで特定されたランダムな効果はモデルから消えなければならないか、ランダムなエフェクトごとにいくつかの値を与えなければなりません。

ここで、DATAステップコードでサポートされていない上記のステートメントの1つは、通常はNLMIXEDコードでも実行する必要があります（！）。特に、BOUNDSステートメントは、NLMIXEDの実行中にあらゆる種類の問題を引き起こす可能性があります。さらに、モデルはほとんど常にBOUNDSステートメントが必要とされないような方法でパラメータ化することができます。負でないことを望む2つのパラメータ（sdn1とsdC）があるので、2つのパラメータ（log_sdn1とlog_sdC）の対数でモデルを書くことができます。 PARMS（またはRETAIN）文の直後に変数sdn1 = exp（log_sdn1）とsdC = exp（log_sdC）を作成してから、sdn1とsdCを尤度の指定で参照することができます。もちろん、BOUNDS文でsdn1とsdCを0にすることができます。 log_sdn1とlog_sdCを指数化すると、0値の（無限に近い）近似しか得られません。

上記で表示した最初の2つのレコードからなるデータセットを作成し、NLMIXEDコードを修正して上記のDATAステップコードに変換しました。下記を参照してください。

data new; 
    input id y1-y5; 
cards; 
1 0.3393279811 2.9297994781 3.3204747713 5.712911709 7.303402636 
2 0.6299412759 0.5226832086 2.6025372538 4.6200632561 6.3933354199 
;;; 

/*proc nlmixed 
     data=new method=GAUSS NOAD qpoints=50 tech=newrap maxit=5000;*/ 
data _null_; 
    set new; 

    array y{*} y1-y5; 
    array Time{*} Time1-Time5; 
    do i=1 to 5; 
    Time{i} = i; 
    end; 

    /* Change PARMS statement to RETAIN statement and express sdn1 and sdC differently */ 
    /*parms b1=0.95 sdn1=0.95 sdC=0.90 fi=0.45;*/ 
    retain  b1 0.95 
      log_sdn1 %sysfunc(log(0.95)) 
      log_sdC %sysfunc(log(0.90)) 
       fi 0.45; 
    /*bounds sdn1 >= 0, sdC >= 0;*/ 
    sdn1 = exp(log_sdn1); 
    sdC = exp(log_sdC); 
    pi=constant('pi'); 

    /* Add bD=0 to data step code only so that there is some value specified for the random effect. */ 
    bD=0; 

    L1 = log(pi) - 
     ((pi*(y[1]-b1*Time[1] + bD))/(sdn1*sqrt(3))) - 
     log(sdn1) - 0.5*log(3) - 
     2*log(1 + exp(-((pi *(y[1]-b1*Time[1] + bD))/(sdn1*sqrt(3))))); 
    L2=1; 
    do j=2 to 5; 
    L2=L2*(log(pi) - ((pi*(y[j]-b1*Time[j]+y[j-1]*fi+ bD))/(sdn1*sqrt(3*(1-fi*fi)))) - 
      log(sdn1) - 0.5*log(3*(1-fi*fi)) - 
      2*log(1+exp(-((pi*(y[j]-b1*Time[j]+y[j-1]*fi+ bD))/(sdn1*sqrt(3*(1-fi*fi)))))) 
     ); 
    end; 
    L=L1*L2; 
    dum=1; 
    /*model dum ~ general(L);*/ 
    /*random bD ~ normal(0,sdC*sdC) subject=ID;*/ 
run; 

最終的なメモをいくつか提供します。私は頻繁に広いデータセットを使用して、対象ごとに少数の観測がある尤度関数を構築してきました。私があなたの尤度関数をどのように構築したかで作った限られた見方から、私は長い形式のデータを再構成することがあなたに利益をもたらすとは思わない。

上記のコードは正常に動作します。しかし、私は2つのレコードだけを使用しました。あなたが除外した他のレコードの可能性関数の問題を特定するかもしれません。

最適化の問題はBOUNDSステートメントの使用に関連している可能性もあります。再パラメータ化されたモデルの使用は、尤度関数が正しく実行されるために必要なすべてのものである可能性があります。