PROC NLMIXED SUMMARY Strengths: - Easy to specify non-linear model - - PDF document

PROC NLMIXED SUMMARY Strengths: - Easy to specify non-linear model - Conditional distribution of Y can be almost anything. Normal, Poisson, Binomial, Gamma, Negative Binomial are built in or you can program your own log- likelihood function - Can estimate non-linear function of parameters with delta method - Easy to generate predicted values with or without EBLUPs Limitations: - Can’t model autocorrelation - Only two levels (Level 1 and Level 2) - Level 2 random effects must be multivariate normal - Tends to be slow for very large problems Using PROC NLMIXED for a normal response with a non-linear model: We use the MIXED.IQ data for an example: Y verbal IQ of post-coma patient SD square root of duration of coma in days DPC time of IQ test in days post-recovery ID identifier for patients 1. First write the model as a formula: Combined (composite) form: − × β = β + + β × + β + ε 0.69 DPC / Y u SD e it h 0 e it i sd i it = + ε EV it it Or more, generally, write a formula for the expected value of Y it − × β = β + + β × + β 0.69 / DPC EV u SD e it h e 0 it i sd i

We can use a multilevel approach: − × β 0.69 / = + β DPC EV B e it h e it i = β + β × + B SD u i 0 i i sd 2. Specify distributions: Level 1: N EV σ 2 | ~ ( , ) Y EV it it it Level 2: τ ~ (0, ) u N 00 i 3. Classify variables 1. Data: Y SD DPC (these are variables in the data set) it i it Response: Y it Predictors: SD DPC i it 2. Basic Parameters: β β β β σ τ 2 e 0 00 sd h 3. Random variables: ε or Y Level 1: it it u Level 2: i 4. Computed variables and parameters: EV B it i 5. OPTIONAL: Variance Model and Parameters: τ positive, we can use a ‘log-parametrization’ in terms e.g. to keep 00 of L, the log of variance. τ = exp( ) L 00 τ will be positive. Note that whatever the value of L may be, 00 2

Changes to the lists 1. Data: same τ , add L : 2. Basic Parameters: drop 00 β β β β σ 2 L e 0 sd h 3. Random: same 4. Computed variables and parameters: τ EV B 00 it i 4. Writing NLMIXED code Preliminary SORT by ID to play it safe: PROC SORT; DATA = IQ; BY ID; RUN; PROC Statement: PROC NLMIXED DATA = IQ QMAX = 300; PARMS Statement: names and initial values of basic parameters: PARMS b0 = 100 bsd = -2.5 be = - 18 s2 = 8 L = 2 bh = 40; Programming statements: Define all computed variables and parameters: B1 = b0 + u; EV = B1 + bsd*SD; + be * exp(-0.6931 * DPC / bh); t2 = exp(L); 3

MODEL Statement: Randomness at Level 1 Distribution of response conditional on random effects MODEL Y ~ normal(EV, s2); RANDOM Statement: Randomness at Level 2: RANDOM u ~ normal(0, t2) SUBJECT = ID; ESTIMATE Statements: Uses Delta method to estimate functions of parameters: ESTIMATE ’IQ at DC=20 DPC = 365’ b0 + bsd * sqrt(20) + be * exp(-0.69*365/bh); ESTIMATE ’Between Sub SD’ exp(L/2); ESTIMATE ’Within Sub SD’ sqrt(s2); ESTIMATE ’Total SD’ sqrt( t2 + s2 ); ESTIMATE ’Reliability’ t2/( t2 + s2 ); PREDICT Statements: PREDICT EV OUT=OUTDS; run; See output in Appendix 4

More than one random effect: Log-Cholesky parametrization of matrix T Suppose we would like to model both the initial deficit and the long-term level as random and depending on the duration of coma. A Log-Cholesky parametrization of T ensures that T is a proper variance matrix. 1, First write the model as a formula (we go directly to the multilevel form): − × β = + 0.69 / DPC EV B B e it h 0 1 it i i = β + β × + B SD u 0 0 0 i 0 sd i i = β + β × + B SD u 1 1 1 i 1 i sd i 2. Specify distributions: Level 1: N EV σ 2 | ~ ( , ) Y EV (same as before) it it it Level 2: ⎛ τ τ ⎞ ⎡ ⎤ ⎡ ⎤ u ⎡ ⎤ 0 ⎜ ⎟ 0 00 01 i ⎢ ⎥ ~ , ⎢ ⎥ ⎢ ⎥ N τ τ ⎜ ⎟ 0 u ⎢ ⎥ ⎜ ⎢ ⎥ ⎟ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ 1 i 01 11 ⎝ ⎠ Log-Cholesky parametrization of T : τ = + 2 2 c c 11 11 10 τ = c c 01 10 00 τ = 2 c 00 00 = exp( ) c L 11 1 = exp( ) c L 00 0 5

3. Classify variables Y SD DPC (same as before) 1. Data: it i it Y Response: it SD DPC Predictors: i it 2. Basic Parameters: Expected value: β β β β β σ 2 L L c 0 1 1 0 01 0 sd 1 sd h Variance: σ 2 L L c 1 0 01 3. Random variables: ε or Level 1: Y it it Level 2: u u 0 1 i i 4. Computed variables and parameters: EV B B it 0 i 1 i τ τ τ c c 00 11 01 00 11 4. Writing NLMIXED code PROC Statement: PROC NLMIXED DATA = IQ QMAX = 300; PARMS Statement: names and initial values of basic parameters: PARMS b0 = 100 b0sd = -2.5 b1 = - 18 b1sd = 0 bh = 40 s2 = 8 L0 = 2 6

L1 = 2 c01 = 0 ; Programming statements: Define all computed variables and parameters: B0i = b0 + b0sd * SD + u0; B1i = b1 + b1sd * SD + u1; EV = B0i + B1i * exp(-0.6931 * DPC / bh); c00 = exp(L0); c11 = exp(L1); t00 = c00**2; t01 = c00*c01; t11 = c11**2 + c01**2; MODEL Statement: Randomness at Level 1 Distribution of response conditional on random effects MODEL Y ~ normal(EV, s2); RANDOM Statement: Randomness at Level 2: RANDOM u0 u1 ~ normal([0, 0], [t00, t01, t11]) SUBJECT = ID; ESTIMATE Statements: Uses Delta method to estimate functions of parameters: ESTIMATE ’IQ at DC=20 DPC = 365’ b0 + b0sd * sqrt(20) + b1 * b1sd * sqrt(20) + exp(-0.69*365/bh); PREDICT Statements: PREDICT EV OUT=OUTDS; run; 7

Using PROC NLMIXED for multilevel logistic regression: In dataset: Y binary outcome X1, X2 predictors ID subject identifier SAS code for logistic regression: PROC NLMIXED DATA=dataset; PARMS b0 = 0 b1 = 0 b2 = 0 L = 0; eta = b0 + b1*X1 + b2*X2 + u; /* linear model */ p = 1 / (1 + exp(- eta )); /* inverse link function */ tau = exp(L); /* variance model */ MODEL Y ~ binary( p ); RANDOM u ~ normal ( 0 , tau ) subject = id; run; 8

Appendix: SAS code for asymptotic regression data iq; set mixed.iq; SD = sqrt(DCOMA); DPC = DAYSPC; Y = VIQ; run; proc contents data = iq; run; proc sort data = iq; by id; run; PROC NLMIXED DATA = IQ QMAX = 300; PARMS b0 = 100 bsd = -2.5 be = - 18 s2 = 8 L = 2 bh = 40; B1 = b0 + u; EV = B1 + bsd*SD + be * exp(-0.6931 * DPC / bh); t2 = exp(L); MODEL Y ~ normal(EV, s2); RANDOM u ~ normal(0, t2) SUBJECT = ID; ESTIMATE 'IQ at DC=20 DPC = 365' b0 + bsd * sqrt(20) + be * exp(-0.69*365/bh); ESTIMATE 'Between Sub SD' exp(L/2); ESTIMATE 'Within Sub SD' sqrt(s2); ESTIMATE 'Total SD' sqrt( t2 + s2 ); ESTIMATE 'Reliability' t2/( t2 + s2 ); PREDICT EV OUT=OUTDS; run; 9

/* model with two random effects */ PROC NLMIXED DATA = IQ QMAX = 300; PARMS b0 = 100 b0sd = -2.5 b1 = - 18 b1sd = 0 bh = 40 s2 = 8 L0 = 2 L1 = 2 c01 = 0 ; B0i = b0 + b0sd * SD + u0; B1i = b1 + b1sd * SD + u1; EV = B0i + B1i * exp(-0.6931 * DPC / bh); c00 = exp(L0); c11 = exp(L1); t00 = c00**2; t01 = c00*c01; t11 = c11**2 + c01**2; MODEL Y ~ normal(EV, s2); RANDOM u0 u1 ~ normal([0, 0], [t00, t01, t11]) SUBJECT = ID; ESTIMATE 'IQ at DC=20 DPC = 365' b0 + b0sd * sqrt(20) + b1 * b1sd * sqrt(20) + exp(-0.69*365/bh); PREDICT EV OUT=OUTDS; run; 10