How to Fit Nonlinear Mixed-Effects Models That Include an - - PowerPoint PPT Presentation
How to Fit Nonlinear Mixed-Effects Models That Include an Intra-Subject (R-side) Covariance Structure Wisconsin-Illinois SAS Users - 2013 November 4, 2013 Edward F. Vonesh (e-vonesh@northwestern.edu) Professor, Department of Preventive
Wisconsin-Illinois SAS Users - 2013 November 4, 2013
Edward F. Vonesh (e-vonesh@northwestern.edu) Professor, Department of Preventive Medicine Feinberg School of Medicine Northwestern University Chicago, IL 60611
1
Supplemental Material Much of the material including programs and data appear in the recently published book: Vonesh, Edward F. 2012. Generalized Linear and Nonlinear Models for Correlated Data: Theory and Applications Using SAS . Cary, NC: SAS Institute Inc. Any examples taken directly from the book will be noted on the slides and likewise any corresponding material including datasets, programs, and SAS macros taken from the book can be found on the author’s web page: http://support.sas.com/publishing/authors/vonesh.html In some cases, the programs presented here have been modified slightly from those given in the book but the essential features remain the same.
2
within-subject (R-side) covariance structure . . . . . . . . . . . . . 18
3
Types of correlated data
each experimental unit is observed on multiple occasions and possibly under different experimental conditions. > Longitudinal data - data in which the underlying metameter (i.e., unit of measurement) for the occasions at which repeated measurements of the response variable are taken is time and interest generally centers on modeling trends over time. > Repeated measurements are more general in that the underlying metameter may be time or it may be a set of experimental conditions (e.g., dose levels of a drug).
4
Types of correlated data (continued)
clusters resulting in within-cluster data that tend to be correlated > Correlation induced by clustering is more often than not accounted for through the specification of a random-effects model in which cluster-specific random effects are used to differentiate within- and between-cluster variability resulting in some form of intra-class correlation. > In some instances, there may be more than one level of clustering resulting in specification of a multi-level random-effects structure. For example, in a multi-center longitudinal trial, we may have a random center effect as well as repeated measurements on subjects within centers.
5
Types of correlated data (continued)
include both a response variable and a location vector associated with that response variable. The location vector describes the position or location at which the measured response was obtained. The proximity of measured responses with one another determines the extent to which they are correlated. > Lattice data, where measurements are linked to discrete regions (e.g., townships, counties, etc.) rather than some continuous coordinate system, are also considered as spatially correlated and are usually obtained from administrative data sources like census data, socio-economical data, and health data.
6
Types of covariates:
measures factors
longitudinal study
7
Examples of repeated measurements:
Clinical studies where the plasma concentration of a drug is measured at several time points for each subject
Econometric studies entailing both cross-sectional and time-series data together in a two-dimensional array.
Clinical studies where each subject receives each treatment on possibly more than one occasion
Pediatric studies examining the growth pattern of children. Agricultural studies examining the growth pattern of plants.
8
Examples of clustered data:
Studies on twins where each pair serves as a natural cluster. Ophthalmology studies where a pair of eyes serves as a cluster.
Studies on members of a litter of animals (toxicology studies). Epidemiology studies of cancer where families serve as clusters.
RCT’s where observations from subjects within a center (cluster) may be correlated (intra-class correlation).
Studies in which different plots of land serve as clusters and measurements within a plot are homogeneous.
9
Examples of spatially correlated data:
Forestry and agronomy studies where sampling from specified (fixed) locations is used to draw inference over an entire region accounting for spatial dependencies.
Studies aimed at describing the incidence and prevalence of a particular disease often rely on spatial correlation to explain regional variations associated with the disease.
Image segmentation studies in the field of medicine are often used to help identify tissue regions that have been stained versus not stained.
10
Linear Models
Generalized Linear Models
Generalized Nonlinear Models
> Normal-theory nonlinear model (NLM)
> Normal-theory nonlinear mixed-effects (NLME) model ◭
11
Nonlinear mixed-effects (NLME) model
The normal-theory nonlinear mixed-effects (NLME) model for correlated responses may be written as yij = f(x′
ij, β, bi) + wij(x′ ij, β, bi)−1/2ǫij, i = 1, . . . , n; j = 1, . . . , pi
where yij is the jth measurement on the ith subject or cluster x′
ij are 1 × u vectors of within- and between-subject covariates
β is an s × 1 unknown vector of regression parameters bi are v × 1 vectors of random-effects ∼ iid N(0, Ψ(θb)) wij(x′
ij, β, bi) are scalar weights that may depend on β and bi
ǫi = (ǫi1, ǫi2, . . . , ǫipi)′ are pi × 1 random error vectors that are mutually independent and distributed as Npi(0, Λi(θw))
12
Nonlinear mixed-effects (NLME) model
We can write this model in matrix form as follows: yi = f(Xi, β, bi) + W i(Xi, β, bi)−1/2ǫi, where yi is the pi × 1 response vector, yi = (yi1 . . . yipi)′ Xi is a pi × u design matrix of within- and between-subject covariates formed from the row vectors x′
ij
ǫi ∼ind Npi(0, Λi(θw)) W i(Xi, β, bi)−1/2 = Diag{wij(x′
ij, β, bi)−1/2}
Let Λi(β, bi, θw) be the pi × pi covariance matrix of (yi|bi) where Λi(β, bi, θw) = W i(Xi, β, bi)−1/2Λi(θw)W i(Xi, β, bi)−1/2 depends on a dw × 1 parameter vector θw and possibly on β and bi.
13
Nonlinear mixed-effects (NLME) model
Using the familar SAS notation from PROC MIXED, we have Λi(θw) is an R-side (within-subject) covariance matrix Ψ(θb) is a G-side (random-effects) covariance matrix
Problem:
While SAS procedures MIXED and GLIMMIX allow for both an R-side and G-side covariance structure that includes directly specifying a within-subject correlation structure, the NLMIXED procedure does not currently support R-side covariance structures. The ǫij are assumed to be conditionally independent given bi
14
Estimation for marginal models (GLIM and GNLM) Generalized Estimating Equations (GEE)
Maximum Likelihood Estimation (MLE)
Estimation for mixed models (GLME and GNLME) Conditional Generalized Estimating Equations (CGEE)
Maximum Likelihood Estimation (MLE)
15
Estimation for mixed models (GLME, GNLME) Conditional GEE (CGEE) based estimation CGEE involves solving a set of first-order CGEE (CGEE1) for (β, b1, . . . , bn) and PL estimating equations for θ = (θb, θw).
This is the approach implemented in the SAS macro %NLINMIX for NLME models and an optional approach in the GLIMMIX procedure for GLME models
16
Estimation for mixed models (GLME, GNLME) Maximum Likelihood (ML) based estimation MLE’s of τ = (β, θ) maximize the integrated log-likelihood L(τ; y) =
n
log{π(yi; β, θ)} =
n
log
where π(yi; β, θw|bi) and π(bi; θb) are the pdf’s of yi|bi and bi.
second-order Laplace approximation to approximate L(τ; y). For the NLME model,
17
include a within-subject covariance structure
when fitting nonlinear mixed-effects (NLME) models
having a conditional within-subject variance-covariance structure Λi(β, bi, θw) = W i(Xi, β, bi)−1/2Λi(θw)W i(Xi, β, bi)−1/2 where W i(Xi, β, bi) is any specified diagonal weight matrix and Λi(θw) is a specified within-subject (R-side) covariance matrix like those found available with PROC MIXED.
18
Within-subject (R-side) covariance structures in MIXED Type Structure: (k, l)th Element θw
Independence
Λi(θw) = σ2Ipi σ2
CS
Λi(θw) = σ2(1 − ρ)Ipi + σ2ρJpi σ2, ρ
AR(1)
Λi(θw) = σ2(ρkl), ρkl = ρ|k−l| σ2, ρ
Spatial (dkl)
Λi(θw) = σ2(ρkl), ρkl = exp{−ρdkl} σ2, ρ
where dkl =Euclidian distance Toeplitz (h+1) Λi(θw) = σ2(ρkl), ρkl =
ρm m ≤ h m > h σ2, ρ m = |k − l|, ρ0 = 1, ρ = (ρ1 . . . ρh)′
Linear (h)
Λi(θw) = θ1Hi1 + . . . + θhHih for θw
known matrices, Hi1, . . . , Hih and unknown parameter vector θw = (θ1 . . . θh)′
19
Fitting nonlinear models in NLMIXED that include a within-subject covariance structure
SOLUTION Consider a class of NLME models where an intra-subject variance-covariance matrix, Λi(θw), can be specified based on a vector of linear random effects. Specifically, let bi = bi1 bi2 be a v × 1 partitioned vector of random effects with v1 × 1 and v2 × 1 component vectors bi1 and bi2, respectively, such that bi1 bi2 ∼ N 0 , Ψ11(θb1) 0′ Ψ22(θb2) . Since bi1 and bi2 are orthogonal to each other, we can form a NLME model which is nonlinear in bi1 but strictly linear in bi2.
20
SOLUTION (continued) Specifically, consider the class of NLME models of the form yi = f(Xi, β, bi1) + W i(Xi, β, bi1)−1/2ǫi where f(Xi, β, bi1) is a general nonlinear function of β and bi1, W i(Xi, β, bi1) = Diag{wij} is a pi × pi diagonal weight matrix with weights wij that may depend on Xi, β and bi1, Zi is a pi × v2 matrix of fixed known constants, ǫi = Zibi2 + ǫ∗
i and ǫ∗ i ∼iid N(0, σ2 wIpi).
Under this NLME model, ǫi is an “intra-subject” error term that is independent of bi1 and which satisfies V ar(ǫi) = Λi(θw) = ZiΨ22(θb2)Z′
i + σ2 wIpi
where θw = (θb2, σ2
w) is a vector of covariance parameters. 21
SOLUTION (continued) By defining select values for Zi and/or σ2
w, one can formulate
Λi(θw) = ZiΨ22(θb2)Z′
i + σ2 wIpi in such a way as to represent
many of the covariance structures available in PROC MIXED. Example 1. Set Zi = 1pi and Ψ22(θb2) = σ2
b2 such that bi2 ∼ N(0, σ2 b2). The
intra-subject covariance matrix Λi(θw) will then satisfy the assumption of compound symmetry in that Λi(θw) = σ2(1 − ρ)Ipi + σ2ρJpi with σ2 = (σ2
w + σ2 b2) and ρ = σ2 b2/(σ2 w + σ2 b2). 22
SOLUTION (continued) Example 2. Assuming one has balanced data with pi = p, one could set Zi = Ip and define Ψ22(θb2) to be any structured p × p matrix. One would then have to set σ2
w to some arbitrarily small value, say
σ2
w = δ ≤ 1E-8 so that Ψ22(θb2) + δIp = Λ(θw) is the desired p × p
intra-subject variance-covariance matrix with θb2 = θw.
mildly unbalanced data with additional programming. The following examples illustrate this approach.
23
Example 1: Bone Mineral Density Data (Lloyd, 1993)
effects of daily calcium supplement in healthy adolescent women.
for the ith subject on the jth occasion.
x1i = 0, if placebo (P) 1, if calcium (C)
24
Bone Mineral Density Data (continued)
TBBMD = β1 + β2 × t
Determine if daily calcium supplementation improves the rate (β2) of bone gain during early adolescence
Fit a cell-means (i.e., full rank) marginal linear model assuming some specified variance-covariance structure and test whether the slopes differ between the two treatment groups.
The dataset is available on the author’s web page but much of the programming presented here is not.
25
Profile plots of TBBMD for adolescent women
26
Bone Mineral Density Data (continued)
yi = Xiβ + ǫi which, for subjects with no missing data, is yi1 yi2 yi3 yi4 yi5 = x1i (1 − x1i) x1iti1 (1 − x1i)ti1 x1i (1 − x1i) x1iti2 (1 − x1i)ti2 x1i (1 − x1i) x1iti3 (1 − x1i)ti3 x1i (1 − x1i) x1iti4 (1 − x1i)ti4 x1i (1 − x1i) x1iti5 (1 − x1i)ti5 β11 β12 β21 β22 + ǫi1 ǫi2 ǫi3 ǫi4 ǫi5 Here, β11 and β21 are the intercept and slope for the calcium group while β12 and β22 are the intercept and slope for the placebo group. We assume ǫi ∼iid N5(0, Σ(θ)).
27
Bone Mineral Density Data (continued)
Assume that for subjects with complete data, ǫi ∼iid N5(0, Σ(θ)) with Σ(θ) = σ2
1
σ1σ2ρ σ1σ3ρ2 σ1σ4ρ3 σ1σ5ρ4 σ2
2
σ2σ3ρ σ2σ4ρ2 σ2σ5ρ3 σ2
3
σ3σ4ρ σ3σ5ρ2 σ2
4
σ4σ5ρ σ2
5
such that we have unequal variances but an AR(1) type correlation across time (heterogeneous AR(1) or TYPE=ARH(1) in MIXED).
1, σ2 2, σ2 3, σ2 4, σ2 5, ρ)′. 28
Bone Mineral Density Data (continued)
autoregressive structure (ARH(1)) using PROC MIXED:
data TBBMD; input (subject group date1 tbbmd1 date2 tbbmd2 date3 tbbmd3 date4 tbbmd4 date5 tbbmd5) (3. @5 $char1. @7 date7. @15 5.3 @21 date7. @29 5.3 @35 date7. @43 5.3 @49 date7. @57 5.3 @63 date7. @71 5.3); if tbbmd1^=. & tbbmd2^=. & tbbmd3^=. & tbbmd4^=. & tbbmd5^=. then data=’Complete ’; else data=’Incomplete’; format date1 date2 date3 date4 date5 date7.;
29
cards;
101 C 01MAY90 0.815 05NOV90 0.875 24APR91 0.911 30OCT91 0.952 29APR92 0.970 102 P 01MAY90 0.813 05NOV90 0.833 15APR91 0.855 21OCT91 0.881 13APR92 0.901 103 P 02MAY90 0.812 05NOV90 0.812 17APR91 0.843 23OCT91 0.855 15APR92 0.895 104 C 09MAY90 0.804 12NOV90 0.847 29APR91 0.885 11NOV91 0.920 19JUN92 0.948 105 C 14MAY90 0.904 10DEC90 0.927 06MAY91 0.952 04NOV91 0.955 27APR92 1.002 106 P 15MAY90 0.831 26NOV90 0.855 20MAY91 0.890 17DEC91 0.908 15JUL92 0.933 107 P 23MAY90 0.777 05DEC90 0.803 15MAY91 0.817 27NOV91 0.809 04MAY92 0.823 108 C 23MAY90 0.792 03DEC90 0.814 109 C 23MAY90 0.821 03DEC90 0.850 06MAY91 0.865 05NOV91 0.879 29APR92 0.908 110 P 30MAY90 0.823 10DEC90 0.827 21MAY91 0.839 25NOV91 0.885 15MAY92 0.923 111 C 30MAY90 0.828 05DEC90 0.873 07AUG91 0.935 12FEB92 0.952 112 P 04JUN90 0.797 17DEC90 0.818 13MAY91 0.817 16DEC91 0.847 18MAY92 0.862 . . . 30
data example2 2 2; set TBBMD; array d date1-date5; array t tbbmd1-tbbmd5; do visit=1 to 5 by 1; date=d[visit]; tbbmd=t[visit]; years=0.5*(visit-1); format date date7.; keep subject group data visit date years tbbmd;
end; run;
31
proc sort data=example2 2 2; by subject visit; run; proc mixed data=example2 2 2 method=ml; class subject group visit; model tbbmd=group group*years/noint solution; repeated visit /subject=subject type=ARH(1); estimate ’slope diff’ group*years 1 -1; run;
32
Covariance Parameter Estimates Cov Standard Z Parm Subject Estimate Error Value Pr Z Var(1) subject 0.003947 0.000527 7.49 < .0001 Var(2) subject 0.004695 0.000630 7.45 < .0001 Var(3) subject 0.004865 0.000660 7.37 < .0001 Var(4) subject 0.005186 0.000712 7.28 < .0001 Var(5) subject 0.004707 0.000646 7.28 < .0001 ARH(1) subject 0.9737 0.003877 251.14 < .0001 Fit Statistics
AIC (smaller is better)
AICC (smaller is better)
BIC (smaller is better)
33
Solution for Fixed Effects Standard Effect group Estimate Error DF t Value Pr > |t| group C 0.8777 0.007979 110 109.99 < .0001 group P 0.8674 0.007836 110 110.70 < .0001 years*group C 0.05304 0.002181 387 24.32 < .0001 years*group P 0.04303 0.002129 387 20.21 < .0001 Estimates Standard Label Estimate Error DF t Value Pr > |t| slope diff 0.01001 0.003047 387 3.29 0.0011 34
How do we fit this model using NLMIXED? First, we need some preliminary matrix notation
2m(m + 1) × 1 column vector from
an m × m matrix A by stacking the lower diagonal elements of A. If A = a11 a12 a13 a21 a22 a23 a31 a32 a33 then V ech(A) = a11 a21 a31 a22 a32 a33 Then we need a SAS macro that performs this operation
35
Bone Mineral Density Data (continued)
first-order autoregressive structure (ARH(1)) using PROC NLMIXED:
/* Macro VECH creates a vech representation of the */ /* intra-cluster covariance matrix defined in the */ /* NLMIXED code for any particular application that*/ /* requires an intra-cluster or within-subject */ /* covariance structure (see Appendix D of Ch 7 */ /* for a more detailed presentation of the macro) */ %macro vech(dim=4, cov=cov, vechcov=vechcov, name=c); %global n vech; %let n=%eval(&dim*(&dim+1)/2);
36
array &vechcov.[&n] &name.1-&name.&n; %let k=0; %do i=1 %to &dim; %do j=1 %to &i; %let k=%eval(&k+1); &vechcov.[&k] = &cov.[&i, &j]; %end; %end; %let temp = &name.1; %do i=2 %to &n; %let temp=&temp,&name.&i; %end; %let vech = &temp; %mend vech;
37
data example2 2 2; set example2 2 2; by subject visit; I1=(visit=1); I2=(visit=2); I3=(visit=3); I4=(visit=4); I5=(visit=5); run; proc nlmixed data=example2 2 2 qpoints=1 tech=newrap; parms b0c=0.8 b0p=.8 b1c=.05 b1p=.05 sig1=0.0004 sig2=0.0004 sig3=0.0004 sig4=0.0004 sig5=0.0004 rho=0.948;
38
eij = I1*e1+I2*e2+I3*e3+I4*e4+I5*e5; predmean = b0c*(group=’C’) + b0p*(group=’P’) + b1c*(group=’C’)*years + b1p*(group=’P’)*years + eij; array cov[5,5]; array sigsq[5] sig1-sig5; do i=1 to 5; do j=1 to 5; dij=1-(i=j); cov[i,j] = sqrt(sigsq[i]*sigsq[j])*(rho**abs(i-j)); end; end; delta =1E-8; %vech(dim=5, cov=cov, name=c); model tbbmd∼normal(predmean, delta );
39
random e1 e2 e3 e4 e5 ∼ normal([ 0, 0, 0, 0, 0 ], [ c1 , c2 , c3 , c4 , c5 , c6 , c7 , c8 , c9 , c10, c11, c12, c13, c14, c15 ]) subject=subject; estimate ’slope diff’ b1c-b1p; run;
40
Fit Statistics
AIC (smaller is better)
AICC (smaller is better)
BIC (smaller is better)
Parameter Estimates Standard Parameter Estimate Error DF t Value Pr > |t| b0c 0.8777 0.008069 107 108.77 < .0001 b0p 0.8674 0.008063 107 107.57 < .0001 b1c 0.05304 0.002183 107 24.30 < .0001 b1p 0.04303 0.002137 107 20.14 < .0001 sig1 0.003946 0.000528 107 7.47 < .0001 41
sig2 0.004695 0.000632 107 7.43 < .0001 sig3 0.004864 0.000662 107 7.34 < .0001 sig4 0.005185 0.000714 107 7.26 < .0001 sig5 0.004706 0.000648 107 7.26 < .0001 rho 0.9737 0.003886 107 250.55 < .0001 Additional Estimates Standard Label Estimate Error DF t Value Pr > |t| slope diff 0.01001 0.003050 107 3.28 0.0014 42
Bone Mineral Density Data (continued) Summary of Results The MIXED and NLMIXED results are nearly identical with 1) minor differences in the fit statistics and parameter estimates 2) minor differences in standard error estimates of θ 3) minor to moderate differences in standard error estimates of β 4) the above differences are due to different computing algorithms MIXED maximizes a profile log-likelihood function with respect to the covariance parameters θ and estimates β via estimated generalized least squares (EGLS) NLMIXED maximizes a joint log-likelihood function for (β, θ) 5) differences in default degree of freedom (DF) values which one can control with DF options available in both procedures
43
Example 2: Orange Tree Data (Draper and Smith, 1981)
y = β1/(1 + exp{−(t − β2)/β3})
44
Profile plots of orange tree growth data
45
Orange Tree Data (continued) Here we fit two random-effects nonlinear logistic growth curve models to the orange tree data (Vonesh, 2012; Chapter 5, §5.4.1). The basic model is given by yij = (β1 + bi1) {1 + exp(−(tij − β2)/β3)} + ǫij (i = 1, . . . , 5; j = 1, . . . , 7) where bi1 ∼ iid N(0, ψ11) are subject-specific random effects Model 1) Conditional independence: Assume ǫi ∼iid N(0, σ2
wIp) independent of bi1
Model 2) Conditional first-order autoregressive, AR(1): Assume ǫi ∼iid N(0, Λi(θw)) independent of bi where θw = (σ2
w, ρ) and the (j, j′)th component of Λi(θw) is defined
by Cov(ǫij, ǫij′|bi1) = σ2
wρ|tij−tij′| so as to account for
unequally spaced observations.
46
Orange Tree Data (continued)
PROC NLMIXED in combination with macro %COVPARMS:
1) help determine the random-effects structure 2) get good starting values
data example4 5 2; input Tree $ Days y; cards; 1 118 30 1 484 58 1 664 87
47
1 1004 115 1 1231 120 1 1372 142 1 1582 145 2 118 33 2 484 69 2 664 111 2 1004 156 2 1231 172 . . . ; proc sort data=example4 5 2 out=example5 4 1; by Tree Days; run;
48
/* Step 1: Fit a fully specified but dummy NLME model */
proc nlmixed data=example5 4 1 qpoints=1; parms b1=175 b2=700 b3=300 sigsq w=10 to 100 by 10; num = b1+bi1; den = 1 + exp(-(Days-(b2+bi2))/(b3+bi3)); predmean = num/den; resid = y - predmean; model y ∼ normal(predmean, sigsq w); random bi1 bi2 bi3 ∼ normal([0,0,0], [0, 0,0, 0,0,0]) subject=Tree;
49
predict predmean out=predout der; id resid; run;
/* Step 2: Run the macro COVPARMS based on NLMIXED output */ %covparms(parms=OLSparms, predout=predout, resid=resid, method=mspl, random=der bi1 der bi2 der bi3, subject=Tree, type=un, covname=psi,
/* Step 3: Final call to macro and final run of NLMIXED */ %covparms(parms=OLSparms, predout=predout, resid=resid, method=mspl, random=der bi1, subject=tree, type=un, covname=psi,
50
proc print data=MLEparms; title ’Starting values from %COVPARMS’; run;
proc nlmixed data=example5 4 1 qpoints=1 ; parms / data=MLEparms; parms num = b1+bi1; den = 1 + exp(-(Days-b2)/b3); predmean = (num/den); model y ∼ normal(predmean, sigsq w); random bi1 ∼ normal(0, psi11) subject=tree; run;
51
1) Select output from GLIMMIX as generated from the initial call to macro %COVPARMS 2) Print of starting values for all parameters based on the second call to %COVPARMS 3) Select output from final call to NLMIXED
1) Select output from GLIMMIX Estimated G Correlation Matrix Effect Row Col1 Col2 Col3 Der bi1 1 1.0000
Der bi2 2
1.0000 Der bi3 3 52
Covariance Parameter Estimates Standard Wald 95% Confidence Cov Parm Subject Estimate Error Bounds UN(1,1) Tree 1324.17 1027.62 442.45 15058 UN(2,1) Tree
1384.03
1355.87 UN(2,2) Tree 17.2882 4702.67 . . UN(3,1) Tree
1127.30
UN(3,2) Tree
2513.66
4111.54 UN(3,3) Tree 4.09E-16 . . . Residual 55.6490 15.7400 34.2273 106.04 2) Starting values from %COVPARMS Obs Parameter Estimate StdError Ratio 1 b1 192.69 19.8035 9.72993 2 b2 728.75 104.67 6.96222 3 b3 353.53 79.9191 4.42360 4 sigsq w 499.43 119.39 4.18333 5 psi11 1007.41 649.99 1.54988 53
3) Select output from NLMIXED Parameters b1 b2 b3 sigsq w psi11 NegLogLike 192.6869 728.7536 353.5302 499.4316 1007.405 149.856998 Iteration History Iter Calls NegLogLike Diff MaxGrad Slope 1 13 133.619673 16.23733 0.113811
2 16 133.528382 0.091291 0.093596
3 18 131.847953 1.680428 0.040079
4 19 131.614493 0.23346 0.016753
5 20 131.58622 0.028272 0.006201
6 22 131.580314 0.005907 0.002407
7 24 131.578286 0.002028 0.00278
8 26 131.572176 0.006109 0.001513
9 28 131.571921 0.000255 0.000379
10 30 131.571888 0.000033 0.000011
11 32 131.571888 4.909E-8 3.837E-6
54
NOTE: GCONV convergence criterion satisfied. Fit Statistics
263.1 AIC (smaller is better) 273.1 AICC (smaller is better) 275.2 BIC (smaller is better) 271.2 Parameter Estimates Standard Parameter Estimate Error DF t Value Pr > t b1 192.05 15.6682 4 12.26 0.0003 b2 727.92 35.2508 4 20.65 < .0001 b3 348.08 27.0810 4 12.85 0.0002 sigsq w 61.5125 15.8824 4 3.87 0.0179 psi11 1003.11 651.53 4 1.54 0.1985 55
Orange Tree Data (continued)
unequally spaced time points using PROC NLMIXED
NLMIXED
/* Step 4. Modify data to run NLMIXED with AR(1) structure */ data example5 4 1; set example5 4 1; by Tree Days;
56
retain index; array Ind[7] I1-I7; if first.Tree then index=0; index+1; do i=1 to 7; Ind[i]=(index=i); end; days1= 118; days2= 484; days3= 664; days4=1004; days5=1231; days6=1372; days7=1582; run; /* Step 5. Add AR(1) correlation to starting values */ data rho; Parameter=’rho’;estimate=.9; run; data MLEparms rho; set MLEparms rho;
57
run; /* Step 6. Fit the NLME model with AR(1) correlation */
proc nlmixed data=example5 4 1 df=4 qpoints=1; parms / data=MLEparms rho; bounds -1<rho< 1; /* The following variable, eij, is a within-tree */ /* error term expressed as a linear combination of */ /* linear random effects that have a first-order */ /* autoregressive structure for unequally spaced */ /* time points as defined by ARRAY statements */ eij = I1*e1+I2*e2+I3*e3+I4*e4+I5*e5+I6*e6+I7*e7; num = b1+bi1; den = 1 + exp(-(Days-b2)/b3);
58
predmean = (num/den) + eij; delta=1e-8; /* The following array statements define the first */ /* order autoregressive covariance matrix assuming */ /* unequally spaced time points within a cluster */ array cov[7,7]; array d[7] days1-days7; do i=1 to 7; do j=1 to 7; cov[i,j] = sigsq w*(rho**abs(d[i]-d[j])); end; end; /* This call to the macro VECH creates a 7x7 lower */ /* diagonal autoregressive covariance matrix with */
59
/* successive elements c1, c2, c3,..., c28 being */ /* the components of the vech representation of */ /* this matrix where c is determined by the macro */ /* variable NAME= option. This is then used in the */ /* specification of a block diagonal covariance */ /* structure where the random effect bi1 is assumed*/ /* orthogonal to (independent of) the linear random*/ /* effects, e1, e2,...,e7 */ %vech(dim=7, cov=cov, name=c); model y ∼ normal(predmean, delta); random bi1 e1 e2 e3 e4 e5 e6 e7 ∼ normal([ 0, 0, 0, 0, 0, 0, 0, 0 ], [ psi11, , c1 ,
60
, c2 , c3 , , c4 , c5 , c6 , , c7 , c8 , c9 , c10, , c11, c12, c13, c14, c15, , c16, c17, c18, c19, c20, c21, , c22, c23, c24, c25, c26, c27, c28 ]) subject=tree; run; proc print data=pe noobs split=’|’; var Parameter Estimate StandardError DF tValue Probt; title ’Final Parameter Estimates’; run; /* Step 7. Compare NLMIXED with SAS macro %NLINMIX */ %nlinmix(data=example5 4 1,
61
procopt=method=ml covtest, model=%str( num = b1+bi1; den = 1 + exp(-(Days-b2)/b3); predv=num/den; ), parms=%str(b1=150 b2=700 b3=350), stmts=%str(class Tree Days; model pseudo y = d b1 d b2 d b3 / noint notest s cl covb; random d bi1 / type=vc subject=Tree; repeated / subject=Tree type=SP(POW)(Days);), expand=zero );
62
1) Select output from NLMIXED with AR(1) correlation 2) Select output from %NLINMIX with AR(1) correlation
1) Select output from NLMIXED with AR(1) correlation Iteration History Iter Calls NegLogLike Diff MaxGrad Slope 1 4 131.575199 0.03834 0.001623
2 6 131.57302 0.002179 0.001535
3 8 131.572169 0.000852 0.001084
4 9 131.571958 0.00021 0.000473
5 10 131.571929 0.000029 0.000095
6 12 131.571929 4.681E-7 0.000021
7 13 131.571927 1.281E-6 0.000126
NOTE: GCONV convergence criterion satisfied. 63
Fit Statistics
263.1 AIC (smaller is better) 275.1 AICC (smaller is better) 278.1 BIC (smaller is better) 272.8 Final Parameter Estimates Standard Parameter Estimate Error DF t Value Pr > t b1 192.05 15.6054 4 12.31 0.0003 b2 727.94 34.1916 4 21.29 < .0001 b3 348.08 26.1946 4 13.29 0.0002 sigsq w 61.4799 16.6643 4 3.69 0.0210 psi11 1007.38 656.43 4 1.53 0.1997 rho 0.8421 512.45 4 0.00 0.9988 64
2) Select output from %NLINMIX with AR(1) correlation Covariance Parameter Estimates Standard Z Cov Parm Subject Estimate Error Value Pr Z d bi1 Tree 1007.37 649.96 1.55 0.0606 SP(POW) Tree 0.8438 448.78 0.00 0.9985 Residual 61.6906 15.9284 3.87 < .0001 Fit Statistics
263.2 AIC (smaller is better) 275.2 Solution for Fixed Effects Standard Effect Estimate Error DF t Value Pr > |t| d b1 192.69 15.7403 28 12.24 < .0001 d b2 728.76 36.0582 28 20.21 < .0001 d b3 353.53 27.3790 28 12.91 < .0001 65
Orange Tree Data (continued) Summary of Results The macro %COVPARMS may be used to 1) identify parameters having associated random effects 2) obtain reasonable starting values for both β and θ. NLMIXED results for Model 1 and Model 2 are nearly identical for the fixed-effects parameters β and variance components (ψ11, σ2
w)
The estimated autoregressive correlation of ρ = 0.8421 would seem to be significant but its standard error estimate of 512.45 (p-value= 0.9988) suggests otherwise 1) the elapsed time (days) between measurements is huge 2) Corr(yi1, yi2) = 0.8421|118−484| = 4.8 × 10−28
66
Orange Tree Data (continued) The NLMIXED and %NLINMIX results for Model 2 are nearly identical with 1) minor differences in the fit statistics and estimates of (β, θ) 2) minor differences in standard error estimates of β and θ Such differences are the result of using different estimation techniques NLMIXED maximizes a joint log-likelihood function for (β, θ) %NLINMIX uses a doubly iterative algorithm which alternates between maximizing a pseudo log-likelihood function with respect to the covariance parameters θ and estimating β via GEE1 (as stipulated by the ’expand=zero’ option) There are also differences in degree of freedom (DF) values which
67
SUMMARY
NLMIXED to obtain MLE’s of the parameters of any GNLME model that includes a within-subject (R-side) variance-covariance structure provided the data are balanced but possibly incomplete (and, in some cases, even when the data are unbalanced).
V ar(yi) = Σi(β, θ) = W i(β)−1/2Σi(θ)W i(β)−1/2 depends on β and Σi(θ) has some specified marginal covariance structure (%NLINMIX yields GEE1 estimates in this case)
68
Key References
Breslow NE and Clayton DG (1993). Approximate inference in generalized linear mixed models. Journal of the American Statistical Association 88:9-25. Draper NR and Smith H (1981). Applied Regression Analysis, 2nd
Godambe VP and Heyde CC (1987). Quasi-likelihood and optimal
Gourieroux C, Monfort A and Trognon A (1984). Pseudo Maximum Likelihood Methods: Theory. Econometrica 52:681-700. Harville DA (1997). Matrix Algebra from a Statistician’s
Henderson HV and Searle SR (1979). Vec and Vech operators for matrices with some uses in Jacobians and multivariate statistics.
69
Canadian Journal of Statistics 7:65-81. Henderson HV and Searle SR (1981). Vec-Permutation matrix, the vec operator and Kronecker products: A review. Linear and Multilinear Algebra 9:271-288. Liang K-Y and Zeger SL (1986). Longitudinal data analysis using generalized linear models. Biometrika 73:13-22. Littell RC, Milliken GA, Stroup WW, Wolfinger RD and Schabenberger O. (2006). SAS for Mixed Models, Second Edition. Cary, NC: SAS Institute Inc. Lloyd T, Andon MB, Rollings N, et. al. (1993). Calcium supplementation and bone mineral density in adolescent girls. Journal of the American Medical Association 270:841-844. Nelson KP, et. al. (2006). Use of the probability integral transformation to fit nonlinear mixed-effects models with nonormal
70
random effects. Journal of Computational and Graphical Statistics 15:39-57. Prentice RL and Zhao LP (1991). Estimating equations for parameters in means and covariances of multivariate discrete and continuous responses. Biometrics 47:825-839. Thall PF and Vail SC (1990). Some covariance models for longitudinal count data with overdispersion. Biometrics 46:657-671. Vonesh EF, Chinchilli VM and Pu K (1996). Goodness-of-fit in generalized nonlinear mixed-effects models. Biometrics 52:572-587. Vonesh EF and Chinchilli VM (1997). Linear and Nonlinear Models for the Analysis of Repeated Measurements. New York: Marcel Dekker, Inc. Vonesh EF, Wang H and Majumdar D (2001). Generalized least squares, Taylor series linearization, and Fisher’s scoring in
71
multivariate nonlinear regression. Journal of the American Statistical Association 96:282-291. Vonesh EF. (2012). Generalized Linear and Nonlinear Models for Correlated Data: Theory and Applications Using SAS . Cary, NC: SAS Institute Inc. Wolfinger RD and Lin X (1997). Two Taylor-series approximation methods for nonlinear mixed models. Computational Statistics and Data Analysis 25:465-490. Zeger SL, Liang KY and Albert PS (1988). Models for longitudinal data: A generalized estimating equation approach. Biometrics 44:1049-1060.
72
Additional example
Example 3: High-flux hemodialyzer data
high-flux hemodialyzers.
ith dialyzer at the jth transmembrane pressure in mmHg (i = 1, . . . , 20 dialyzers; 10 dialyzers per each of two blood flow rates; j = 1, . . . , 7 measurements per dialyzer).
for the ith dialyzer on the jth occasion with target values of TMP set at 25, 50, 100, 150, 200, 250, 300 mmHg.
73
denoted by ai = 0, if Qb = 200 ml/min 1, if Qb = 300 ml/min (SAS variable: QB 300)
Structurally characterize the ultrafiltration profile of the typical dialyzer as a function of transmembrane pressure and blood flow rate and select a suitable covariance structure. As such, this provides another example for how one can fit nonlinear models in NLMIXED that include a within-subject (R-side) covariance structure.
74
High-flux hemodialyzer data (continued)
model to the data. Preliminary analyses suggest the following NLME model may be used. yij = βi1(1 − exp{−βi2(xij − βi3)}) + ǫij, (i = 1, . . . , 20; j = 1, . . . , 7) βi1 = (β11 + β12ai) + (β13 + β14ai)xij + bi1 βi2 = β21 + β22ai + bi2 βi3 = β31 + β32ai. Here the asymptote parameter, βi1, is a linear function of both blood flow rate (ai) and transmembrane pressure (xij). Initial analysis assumed ǫij ∼iidN(0, σ2
w), i.e., homogeneous
within-dialyzer variation.
75
High-flux hemodialyzer data (continued)
model, it is not clear whether the observed heterogeneity in the data can be adequately explained through specification of the random effects bi1 and bi2 alone or whether a better fit to the data is possible under some alternative covariance structure. To that end, several covariance structures were fit to the data including some of the marginal covariance structures evaluated by Littell et.
rather than fit a quartic polynomial growth curve model assuming different covariance structures as was done by Littell et. al. (2006), Vonesh (2012) fit both mixed-effects and marginal versions of the above linear asymptotic growth curve model assuming different covariance structures.
76
High-flux hemodialyzer data (continued)
covariance structures, Λi(β, bi, θw). Model Λi(β, bi, θw) Intra-Subject Covariance Structure 1 VC constant variance component 2 VC-POM constant variance, power of mean 3 CS-POM compound symmetric, power of mean 4 AR(1)-POM first-order autoregressive, power of mean 5 CSH heterogeneous compound symmetric 6 ARH(1) heterogeneous first-order autoregressive
Key: Models 1-4 include bi1, bi2 as random effects; Models 5-6 include
bi1 as the only random effect.
77
High-flux hemodialyzer data (continued)
Model Σi(β, θ) Marginal Covariance Structure 1 UN unstructured variance-covariance matrix 2 CS compound symmetric 3 CSH heterogeneous compound symmetric 4 CS-POM compound symmetric, power of mean 5 AR(1) first-order autoregressive 6 ARH(1) heterogeneous first-order autoregressive 7 AR(1)-POM first-order autoregressive, power of mean
78
High-flux hemodialyzer data (continued)
together with code similar to that shown for the Orange Tree Data (see slide 44) that allows one to fit various R-side covariance structures.
the following slide. Smaller values of the likelihood-based information criterion (AIC, AICC, BIC) are indicative of a better
Note: Several attempts to fit model 4 under the class of NLME models failed and so this model is not summarized in the following
79
Output 30: Rank order of model goodness-of-fit statistics
Information Criterion Covariance AIC AICC BIC Model Structure Value Rank Value Rank Value Rank 1 NLME 1 RE(bi1,bi2)|VC 32.27 7 34.72 6 44.21 6 NLME 2 RE(bi1,bi2)|VC-POM 24.45 5 27.34 4 37.39 4 NLME 3 RE(bi1,bi2)|CS-POM 24.70 6 28.06 5 38.64 5 NLME 5 RE(bi1)|CSH 42.47 9 47.49 9 59.40 9 NLME 6 RE(bi1)|ARH(1) 6.38 2 11.40 2 23.31 3 2 GNLM 1 UN 15.35 4 41.22 7 51.20 7 GNLM 2 CS 77.58 12 79.28 12 87.54 12 GNLM 3 CSH 54.13 11 58.55 11 70.06 11 GNLM 4 CS-POM 52.57 10 54.63 10 63.52 10 GNLM 5 AR(1) 41.37 8 43.07 8 51.33 8 GNLM 6 ARH(1) 4.91 1 9.34 1 20.84 1 GNLM 7 AR(1)-POM 11.09 3 13.15 3 22.04 2 80
High-flux hemodialyzer data (continued) Summary of Results
GNLM Model 6 with a heterogeneous first-order autoregressive structure, ARH(1), provided the best fit to the data followed closely by either the NLME Model 6 with random asymptote bi1 and intra-subject ARH(1) structure and the GNLM Model 7 with a first-order autoregressive structure coupled with a power of the mean variance (AR(1)-POM).
similar models considered by Littell et. al. (2006) when they fit LM’s and LME models (polynomial growth curves) to this data (SAS for Mixed Models, Chapter 9, §9.5, pp. 374-392).
81
82