RE: Fitting Data Below the Quantification Limit
From: Lewis B. Sheiner - lbs@lewisbsheiner.net
Subject: RE: [NMusers] Fitting Data Below the Quantification Limit
Date: 1/10/2004 5:21 PM
All -
As I understand it, NM6 should have built-in code to do the necessary
integrations in a conveneint way. In the mean-time, here (below and attached)
is some code that Russ Wada, Bill Poland, and I came up with that works at least
on a very simple test problem. That's all the testing I've done. The code assumes Y|eta
is normally distributed with expected value = MN and standard deviation = SD, so it will
not work for time-to-event Y, logisticY, etc.
LBS.
==========================================================
$ERROR [or $PRED]
MN= ... ; Your model for E(DV); often just F
SD= ... ; Your model for SD(DV)
QL= ... ; The QL -- A fixed value
X=(QL-MN)/SD
ABSX=X
IF (X.LT.0) ABSX=-X
; compute -2LL contrib
l2pi=1.837877 ; log 2*pi
IF (DV .GT. QL) THEN ; >QL case: L=N(MN,SD2)
Y = l2pi + 2*LOG(SD)+((DV-MN)/SD)**2
ENDIF
IF (DV .LE. QL) THEN; >QL case: L=N(MN,SD2), need CUM NL
; Compute NN = Std CUM NL (STDNCDF) at X -- from Abramovitz & Stegun,
; Handbook of Mathematical Functions
b1=0.31938153
b2=-0.356563782
b3=1.781477937
b4=-1.821255978
b5=1.330274429
pp=0.2316419
c2=0.3989423
tt=1/(1+ABSX*pp)
b=c2*EXP(-1/2*ABSX**2)
n1=((((b5*tt+b4)*tt+b3)*tt+b2)*tt+b1)*tt
n2=1-b*n1
ENDIF
IF(DV .LE. QL.AND.X.LT.0) THEN ; belongs in above 'IF' but can't nest
NN=1-n2
ELSE
NN = n2
ENDIF
IF (DV.LT.QL.AND.ABSX.LT.6) THEN
Y = -2*LOG(NN)
ENDIF
IF(DV .LE. QL .AND. X .GE. 6.0) THEN ; X > +6SD, STDNCDF~ 1
Y = 0
ENDIF
IF(DV .LE. QL .AND. X .LE. -6.0) THEN ; X < -6SD, trunc to -2*log(1-STDNCDF(6))
Y = -37.8379 ;
ENDIF
$EST METHOD=COND LAPLACE -2LL
$THETA ...
$OMEGA ...
====================================================================
_______________________________________________________