Stat 504, Lecture 28 1✬✫✩✪Non-LinearMixed-Effect ModelsLast time, we began to look at linear mixed-effectmodels for longitudinal data. The model for a singlesubject i with niresponse occasions isyi= Xiβ + Zibi+ i,i=1,...,mwhereyi=(yi1,yi2,...,yi,ni)Tbi∼ Nq(0,ψ)i∼ Nni(0,σ2Vi)β = fixed effects for populationbi= random effects for subject iψ = between-subject covariance matrixσ2Vi= within-subject covariance matrixThe times of measurement ti1,ti2,...,tiniaretypically incorporated into the columns of Xiandpossibly Zi.Stat 504, Lecture 28 2✬✫✩✪We discussed how to fit this type of model in PROCMIXED. This approach is fully parametric, and theresulting ML estimates are appropriate for a responsethat is continuous and approximately normallydistributed.What can we do for longitudinal responses that areclearly non-normal and discrete? The naturalextension is to allow the responses to come from anexponential-family distribution, and to connect themean of this distribution to subject i’s linearpredictor Xiβ + Zbiby a link function. The resultingclass is called:• nonlinear mixed-effects models• generalized linear mixed models• generalized linear models with random effectsLike GLIM’s, these models are regarded as linear ornonlinear, depending on your point of view.Stat 504, Lecture 28 3✬✫✩✪Example: binary responsesResponse is yi=(yi1,yi2,...,yi,ni)T, whereyij= 1 if subject i displayed the traitof interest at occasion j, 0 otherwise∼ Bernoulli(πij)ηij= logit(πij)=log[πij/ (1 − πij)]ηi=(ηi1,ηi2,...,ηi,ni)T= Xiβ + Zibibi∼ N(0,ψ)With a binary-response model, it’s difficult to haveany random effects beyond a random interceptZi=(1, 1,...,1)Tif lots of subjects haveyi=(0, 0,...,0)Tor (1, 1,...,1)T,because slopes for these subjects are poorly estimated.Stat 504, Lecture 28 4✬✫✩✪The likelihood functionL =iP (yi|bi) P (bi) dbicannot be computed exactly; it can only beapproximated by numerical techniques• Gauss-Hermite quadrature• adaptive quadrature• Laplace expansionsAlgorithms for model fitting are considerably morecomplicated than for the normal linear mixed modelStat 504, Lecture 28 5✬✫✩✪Model fitting• Early programs used “penalized quasilikelihood”(PQL) (Breslow and Clayton, 1993).– PQL estimates can be badly biased for thebinary model– even with bias corrections (Lin and Breslow,1996), biases are still bad when variance ofrandom intercepts is large (say, 0.5 or greater)– PQL is no longer recommended• Later programs use true ML, which is muchbetter– HLM– MIXOR (designed for ordinal responses, butalso works for binary)– MLWin– PROC NLMIXED• Bayesian inference is also possible by Markovchain Monte Carlo (MCMC)– MLwin– WinBUGSStat 504, Lecture 28 6✬✫✩✪Modeling issues1. Elements of β estimate subject-specific (SS) butnot population-average (PA), effects. SS and PAhave different interpretations and are appropriatein different circumstances (Zeger, Liang andAlbert, 1988; Neuhaus, Kalbfleisch and Hauck,1991; Hu et al., 1998). The difference between SSand PA grows as between-subjects variationincreases.2. Assuming an underlying normal distribution forthe random effects is often unrealistic; there maybe a large group ofyi=(0, 0,...,0)Tsubjects for which πij≈ 0 at all occasions(Carlin, Wolfe, Brown and Gelman, 2001)Missing dataAs long as you use true ML or Bayesian techniques(not PQL), the implications of missing responses areno different than for normal linear mixed models; theprocedures work as long as the missing values are“missing at random.”Stat 504, Lecture 28 7✬✫✩✪Example. Here is a hypothetical dataset containing5 repeated binary measurements for a sample of 25subjects.WeekSubject012341110112111113000004000005100016011107010118101119111111000011110111112100111300100140000115111111611111170001118110111901101200001021011112200111230100124010002501001Stat 504, Lecture 28 8✬✫✩✪Plot of the sample proportion of successes by week:012340.0 0.2 0.4 0.6 0.8 1.0weekproportion•••••I created these data by supposing that the responsefor subject i at week j is Bernoulli withlog πij1 − πij= αi+ γweekjwhere γ = .5 andαi∼ N(−1, 2).That is, each subject’s logit-probabilities of successfollow a linear trajectory with a random intercept andfixed slope of 0.5; the subjects’ intercepts arerandomly distributed about −1.0 with variance 2.0Stat 504, Lecture 28 9✬✫✩✪In the notation of a nonlinear mixed-effects model,the response vector for subject i,yi=(yi1,yi2,yi3,yi4,yi5)Thas meanµi=(µi1,µi2,µi3,µi4,µi5)T(in this case, the µij’s are actually πij’s). The logit ofthis vector isηi= logit(µi)=Xiβ + Zibi,whereXi=1011121314,Zi=11111,β=−1.00.5,andbi∼ N(0, 2).Stat 504, Lecture 28 10✬✫✩✪The interpretation of the elements of β is a bit subtle.The intercept β0is an average log-odds of success atweek 0 (“average” means averaged across subjects).But exp(β0)isnot the average odds of success atweek 0 in the population, because the average of thelog is not the log of the average. As long as thecoefficients vary from one subject to another, theaverage of the log won’t be the log of the average.For this reason, the elements of β are calledsubject-specific (SS) rather than population-average(PA) effects.Let’s try to fit this model in SAS using PROCNLMIXED. First, we arrange the data in a file calledfake.dat.101111120131141201211221231241300-lines omitted -25 0 025 1 125 2 025 3 025 4 1Stat 504, Lecture 28 11✬✫✩✪Here’s the SAS code.options nocenter nodate nonumber linesize=72;data fake;infile "d:\jls\stat504\lectures\fake.dat";input id week y;run;proc nlmixed data=fake;parms beta0=-1.0 beta1=0.5 varb0=2.0;eta = beta0 + beta1*week + b0;mu = exp(eta) / ( 1 + exp(eta) );model y ~ binary(mu);random b0 ~ normal(0,varb0) subject=id;run;The params statement gives initial values for theparameters to start the interative fitting procedure.NLMIXED requires good starting values! If you don’tgive a params statement, all the starting values areset to 1.0, which can be very bad. In this example, Iused the true values for the parameters as the startingvalues.The next four lines after params define the model.The types of models that can be fit include