topic Re: Nonlinear Mixed Effect Model Analysis in Statistical Procedures

Nonlinear Mixed Effect Model Analysis

Sandhu — Thu, 27 Dec 2012 12:07:22 GMT

Hello.. I am doing analysis of repeated measure data obtained from a clinical trial. There are total 12 time points (one baseline and 11 follow-up) and time is coded in Hours with baseline time coded as "0". Outcomes of trial are three, one continuous (pain assessment using 100 mm VAS scale) and two binary (both binary outcomes coded as 0 and 1). Two predictors (covariates) are sex (male coded as 0 and female coded as 1) and age groups (pre-adolescent coded as 0 and post-adolescent coded as 1). Therefore from statistical analysis point of view, it is a 2x2 factorial design. Objective of my trial are to evaluate the main effect and interaction effect of age and sex on pain (VAS score). The mean profile plot shows a nonlinear quadratic trend of VAS score. I am familiar with PROC MIXED AND PROC GLIMMIX procedure and one possibility is to use polynomial random coefficient model analysis using these two aforementioned procedures. But my objective is not to investigate the trend of pain but to evaluate the overall effect of age, sex and age sex interaction on VAS score. Therefore, i was just thinking about the use of PROC NLMIXED in this case. But honestly speaking, i find it very difficult to use NLMIXED procedure and not able to understand syntax and its application. for example, how to get initial PARMS values etc etc. PROC GLIMMIX is very versatile procedure and i am using it for my two binary outcomes. Is there any possibility that i can use PROC GLIMMIX for nonlinear longitudinal data analysis instead of using PROC NLMIXED?

Regards

Re: Nonlinear Mixed Effect Model Analysis

SteveDenham — Thu, 27 Dec 2012 12:48:37 GMT

When I think of nonlinear, I assume it to mean nonlinear in the parameters, such that a polynomial is not what I think of as nonlinear. Given that, it makes perfect sense to me to use GLIMMIX to explore this. If you have access to SAS/STAT 12.1, check out the EFFECT statement. It will enable you to examine polynomial effects specifically (as well as spline effects, which may be more appropriate).

Steve Denham

Re: Nonlinear Mixed Effect Model Analysis

Sandhu — Thu, 27 Dec 2012 14:24:10 GMT

Dear Steve Denham,

Thanks a lot for replying. I am using SAS 9.2 and do not have access to SAS/STAT 12.1. Following is the code that i can think of fitting a second order (quadratic) polynomial trend by using GLIMMIX

proc glimmix data=data ORDER=Data nobound;

class id sex age;

model vas = time time*time

sex age sex*age

sex*time age*time sex*age*time

sex*time*time age*time*time sex*age*time*time

/ ddfm=kr s cl ;

random int time time*time / subject=id type=un v;

lsmeans sex / pdiff adjust=simulate(acc=.0005 seed=121211 report) adjdfe=row cl;

lsmeans age / pdiff adjust=simulate(acc=.0005 seed=121211 report) adjdfe=row cl;

lsmeans sex*age / pdiff adjust=simulate(acc=.0005 seed=121211 report) adjdfe=row cl;

run;

is this code appropriate for analysis of quadratic trend? should i use random _residual_ statement also to model covariance structure e.g.

random _residual / sub=id type=sp(pow)(xtime). i am thinking of spatial power because my time points are not equally spaced. i can use xtime as class variable and time as continuous variable. So code will be like this:

proc glimmix data=data ORDER=Data nobound;

class id sex age xtime;

model vas = time time*time

sex age sex*age

sex*time age*time sex*age*time

sex*time*time age*time*time sex*age*time*time

/ ddfm=kr s cl ;

random int time time*time / subject=id type=un v;

random _residual / sub=id type=sp(pow)(xtime);

lsmeans sex / pdiff adjust=simulate(acc=.0005 seed=121211 report) adjdfe=row cl;

lsmeans age / pdiff adjust=simulate(acc=.0005 seed=121211 report) adjdfe=row cl;

lsmeans sex*age / pdiff adjust=simulate(acc=.0005 seed=121211 report) adjdfe=row cl;

run;

is this code for analysis and multiple comparison by using lsmeans statement appropriate? or should i use ESTIMATE statement for doing multiple comparisons.

Thanks a lot for the help

Regards

Sandhu

Re: Nonlinear Mixed Effect Model Analysis

SteveDenham — Thu, 27 Dec 2012 14:37:12 GMT

This looks exactly like what you will need. I will be curious to find out if the random _residual_ statement will work. I think it should, but I worry some that if you have a large number of xtime levels, you may have convergence problems.

Good luck, and please let the community know if this approach worked.

Steve Denham

Re: Nonlinear Mixed Effect Model Analysis

Sandhu — Fri, 28 Dec 2012 03:46:43 GMT

Dear Steve Denham,

I tried this code with random _residual_ statement for my data which comprised of 280 subject and each subject had 12 repeated measurements. your prediction was correct. It did not work. The log was :

NOTE: Convergence criterion (ABSGCONV=0.00001) satisfied.

WARNING: Optimization stopped because of infinite objective function.

NOTE: PROCEDURE GLIMMIX used (Total process time):

real time 1.01 seconds

cpu time 0.96 seconds

It converged but "Optimization stopped because of infinite objective function." this happened with both Dual Quasi-Newton and Newton-Raphson with Ridging optimization techniques. If i am correct, the Newton-Raphson with Ridging technique makes PROC GLIMMIX almost similar to PROC MIXED. So this mean random _residual_ in PROC GLIMMIX or Repeated statement Repeated xtime/ subject=id type=sp (pow)(xtime) in PROC MIXED will not work? if so, then how to model covariance structure.

Regards

Sandhu

Re: Nonlinear Mixed Effect Model Analysis

SteveDenham — Fri, 28 Dec 2012 19:49:19 GMT

If the model converged but the optimization stopped, I would look at the ddfm=kr option on the model statement as the most likely suspect, as this is a secondary process applied to the converged model. I am pretty sure that the unequal spacing leads to some problems with the full Kenward-Rogers procedure under unbalanced situations. You might be able to work around this with ddfm=kr(firstorder). If your data are balanced across sex and age, then you could try ddfm=satterthwaite. And if there are still problems, ddfm=bw, where bw is between-within.

Steve Denham

Re: Nonlinear Mixed Effect Model Analysis

Sandhu — Sat, 29 Dec 2012 07:05:35 GMT

Dear Steve Denham,

There are couple of things which i feel are responsible for these problems.

When i was running the original quadratic code, i was getting G matrix not positive warning as follow:

proc glimmix data=data ORDER=Data ;

class id sex age ;

model vas = time time*time

sex age sex*age

sex*time age*time sex*age*time

sex*time*time age*time*time sex*age*time*time / ddfm=kr s cl ;

random int time time*time/ subject=id type=un v;

run;

NOTE: Convergence criterion (GCONV=1E-8) satisfied.

NOTE: At least one element of the gradient is greater than 1e-3.

NOTE: Estimated G matrix is not positive definite.

NOTE: PROCEDURE GLIMMIX used (Total process time):

real time 22.37 seconds

cpu time 22.12 seconds

Therefore, i used nobound option in the procedure as recommended (by Little 2006) for PROC MIXED procedure:

proc glimmix data=data ORDER=Data nobound ;

class id sex age ;

model vas = time time*time

sex age sex*age

sex*time age*time sex*age*time

sex*time*time age*time*time sex*age*time*time / ddfm=kr s cl ;

random int time time*time/ subject=id type=un v;

run;

It worked fine for quadrartic code even after adding nobound.

Now when i added random _residual_ statement in this above code along with nobound option, the problem of "Optimization stopped because of infinite objective function." as i mentioned in the previous post. This problem persisted even with ddfm=kr(firstorder), ddfm=satterthwaite and ddfm=bw .

But when i removed nobound, it worked even with ddfm=kr.

proc glimmix data=data2 ORDER=Data ;

class id sex age xtime;

model vas = time time*time

sex age sex*age

sex*time age*time sex*age*time

sex*time*time age*time*time sex*age*time*time / ddfm=kr s cl ;

random int time time*time/ subject=id type=un v;

random _residual_ / sub=id type=sp(pow)(xtime);

run;

Log is:

NOTE: Convergence criterion (GCONV=1E-8) satisfied.

NOTE: At least one element of the gradient is greater than 1e-3.

NOTE: Estimated G matrix is not positive definite.

NOTE: PROCEDURE GLIMMIX used (Total process time):

real time 23.30 seconds

cpu time 23.10 seconds

Now here is my dilemma : should i use nobound option and no random _residual_ statement OR should i use random _residual_ statement and no nobound option. IF i am not using nobound option, then how to justify the effect of Estimated G matrix is not positive definite.

These above codes were with default Optimization Technique i.e. Dual Quasi-Newton.

Now part two of the problem is changing Optimization Technique from Dual Quasi-Newton to Newton-Raphson with Ridging with nloptions tech= nrridge.

when i run same quadratic code (without random _residual_ statement) along with nobound option as :

proc glimmix data=data2 ORDER=Data nobound;

class id sex age xtime;

model vas = time time*time

sex age sex*age

sex*time age*time sex*age*time

sex*time*time age*time*time sex*age*time*time / ddfm=kr s cl ;

random int time time*time/ subject=id type=un v;

nloptions tech= nrridge; *Newton-Raphson with Ridging technique*;

run;

there is no difference in results for from previous Optimization Technique i.e Dual Quasi-Newton. The Fit criteria, Type III test of fixed effect are exactly same, so as the :

NOTE: Convergence criterion (GCONV=1E-8) satisfied.

NOTE: Estimated G matrix is not positive definite.

NOTE: PROCEDURE GLIMMIX used (Total process time):

real time 3.75 seconds

cpu time 3.71 seconds

This Newton-Raphson with Ridging technique also did not work with nobound option as it happened with Dual Quasi-Newton after i added random _residual_ statement.

BUT The problem occurs when i use nloptions tech= nrridge (i.e Newton-Raphson with Ridging technique) along with random _ residual_ statement but without nobound option.

proc glimmix data=data ORDER=Data ;

class id sex age xtime;

model vas = time time*time

sex age sex*age

sex*time age*time sex*age*time

sex*time*time age*time*time sex*age*time*time / ddfm=kr s cl ;

random int time time*time/ subject=id type=un v;

random _residual_ / sub=id type=sp(pow)(xtime);

nloptions tech= nrridge;

run;

now log is:

NOTE: Convergence criterion (XCONV=0) satisfied.

NOTE: At least one element of the gradient is greater than 1e-3.

NOTE: PROCEDURE GLIMMIX used (Total process time):

real time 10.43 seconds

cpu time 10.34 seconds

IN CONTRAST to log from same statement but with Dual Quasi-Newton technique, the above log (with Newton-Raphson with Ridging technique) shows No NOTE of Estimated G matrix is not positive definite and Convergence criterion changes from (GCONV=1E-8) to (XCONV=0). Another important differences are in Fit Statistics and Type III tests of fixed effects

Fit Statistics (from Dual Quasi-Newton technique)

-2 Res Log Likelihood 24037.76

AIC (smaller is better) 24053.76

AICC (smaller is better) 24053.80

BIC (smaller is better) 24082.84

CAIC (smaller is better) 24090.84

HQIC (smaller is better) 24065.42

Generalized Chi-Square 424340.5

Gener. Chi-Square / DF 126.74

Fit Statistics ( Newton-Raphson with Ridging technique)

-2 Res Log Likelihood 28283.79

AIC (smaller is better) 28299.79

AICC (smaller is better) 28299.83

BIC (smaller is better) 28328.87

CAIC (smaller is better) 28336.87

HQIC (smaller is better) 28311.45

Generalized Chi-Square 939474.5

Gener. Chi-Square / DF 280.61

Type III Tests of Fixed Effects (from Dual Quasi-Newton technique)

Num Den

Effect DF DF F Value Pr > F

TIME 1 1048 90.45 <.0001

TIME*TIME 1 2540 0.84 0.3598

SEX 1 131.9 4.35 0.0389

AGE 1 131.9 1.97 0.1626

SEX*AGE 1 131.9 5.22 0.0240

TIME*SEX 1 1048 0.03 0.8587

TIME*AGE 1 1048 0.00 0.9442

TIME*SEX*AGE 1 1048 0.09 0.7703

TIME*TIME*SEX 1 2540 0.10 0.7469

TIME*TIME*AGE 1 2540 0.01 0.9334

TIME*TIME*SEX*AGE 1 2540 0.00 0.9970

Type III Tests of Fixed Effects ( Newton-Raphson with Ridging technique)

Num Den

Effect DF DF F Value Pr > F

TIME 1 537 13.28 0.0003

TIME*TIME 1 3348 0.76 0.3819

SEX 1 135.9 0.10 0.7577

AGE 1 135.9 0.05 0.8315

SEX*AGE 1 135.9 0.12 0.7270

TIME*SEX 1 537 0.14 0.7074

TIME*AGE 1 537 0.03 0.8522

TIME*SEX*AGE 1 537 0.03 0.8736

TIME*TIME*SEX 1 3348 0.00 0.9580

TIME*TIME*AGE 1 3348 0.00 0.9760

TIME*TIME*SEX*AGE 1 3348 0.00 0.9762

WHY is such as dramatic difference between these two technique when incorporate covariance structure (i.e. random _residual_ statement) with spatial power structure. (I am not sure same effect will be observed if we use ar(1) or any other covariance strucute, BUT definitely it happened with sp(pow)(xtime) covariance structure) AND what is this criterion difference from (GCONV=1E-8) to (XCONV=0).

According to LITTLE 2006, Newton-Raphson with Ridging technique makes PROC GLIMMIX almost identical to PROC MIXED.

I really feel that i am asking very exhausting questions. Thanks for your patience and all the help.

Regards

Sandhu

Re: Nonlinear Mixed Effect Model Analysis

SteveDenham — Mon, 31 Dec 2012 13:21:36 GMT

Without extensive investigation I can't guarantee that this is what is occurring, but I think the NRRIDG technique is finding a saddlepoint in the response surface. To investigate this, you might save the covariance parameter estimates from the QUANEW and use them as starting values.

But what is very intriguing is the large change in the denominator degrees of freedom. This is undoubtedly due to the Kenward-Rogers adjustment, and QUANEW resulting in at least one of the G-side covariance parameters being set to zero. Which parameter has a zero estimate? Can it be removed, and have the model still make sense? If it is the intercept term, consider that your R-side model may be exhausting all of the variation, leaving none to be estimated by the intercept. If this is the case, I think you can safely remove "int" from the first random statement. If not, could you please post all of the variance parameters from both techniques?

Steve Denham

Re: Nonlinear Mixed Effect Model Analysis

Sandhu — Thu, 03 Jan 2013 13:39:14 GMT

Dear Steve Denham,

I can not figure it out that what is causing such a large change with use of these two techniques. I think you are the best person to understand it. Therefore, i am sending code, model information, dimensions, Optimization Information, LOG etc. Interestingly, When i remove intercept, the NRRIDG technique does not work and following error is reported ERROR: "NRRIDG Optimization cannot be completed. Optimization routine cannot improve the function value" whereas in same condition i.e. without intercept, QUANEW technique works and it converges.

A) CODE (NRRIDG with random intercept)

proc glimmix data=data2 ORDER=Data ;

class id sex age xtime;

model vas = time time*time

sex age sex*age

sex*time age*time sex*age*time

sex*time*time age*time*time sex*age*time*time / ddfm=kr s cl ;

random int time time*time/ subject=id type=un v;

random _residual_ / sub=id type=sp(pow)(xtime);

nloptions tech= NRRIDG;

run;

Model Information

Data Set WORK.DATA2

Response Variable VAS

Response Distribution Gaussian

Link Function Identity

Variance Function Default

Variance Matrix Blocked By ID

Estimation Technique Restricted Maximum Likelihood

Degrees of Freedom Method Kenward-Roger

Fixed Effects SE Adjustment Kenward-Roger

Dimensions

G-side Cov. Parameters 6

R-side Cov. Parameters 2

Columns in X 27

Columns in Z per Subject 3

Subjects (Blocks in V) 280

Max Obs per Subject 12

Optimization Information

Optimization Technique Newton-Raphson with Ridging

Parameters in Optimization 7

Lower Boundaries 4

Upper Boundaries 1

Fixed Effects Profiled

Residual Variance Profiled

Starting From Data

Iteration History

Objective Max

Iteration Restarts Evaluations Function Change Gradient

0 0 4 30662.837656 . 8506.067

1 0 2 29840.553145 822.28451046 4711.142

2 0 2 29369.643328 470.90981779 3582.741

3 0 2 28916.249747 453.39358027 3044.089

4 0 3 28783.541782 132.70796514 2944.675

5 0 2 28280.91676 502.62502243 8416.694

6 1 5 28751.288165 -470.3714047 7813.525

Convergence criterion (XCONV=0) satisfied.

Covariance Parameter Estimates

Standard

Cov Parm Subject Estimate Error

UN(1,1) ID 52.0689 .

UN(2,1) ID 1.0951 .

UN(2,2) ID 25.7652 .

UN(3,1) ID -2.3731 .

UN(3,2) ID -10.2573 .

UN(3,3) ID 12.3989 .

SP(POW) ID 0.5592 0.03979

Residual 280.30 32.7243

LOG

NOTE: Convergence criterion (XCONV=0) satisfied.

NOTE: At least one element of the gradient is greater than 1e-3.

NOTE: PROCEDURE GLIMMIX used (Total process time):

real time 10.49 seconds

cpu time 10.43 seconds

B) CODE (NRRIDG without random intercept)

proc glimmix data=data2 ORDER=Data ;

class id sex age xtime;

model vas = time time*time

sex age sex*age

sex*time age*time sex*age*time

sex*time*time age*time*time sex*age*time*time / ddfm=kr s cl ;

random time time*time/ subject=id type=un v;

random _residual_ / sub=id type=sp(pow)(xtime);

nloptions tech= NRRIDG;

run;

Model information remains same

Dimensions

G-side Cov. Parameters 3

R-side Cov. Parameters 2

Columns in X 27

Columns in Z per Subject 2

Subjects (Blocks in V) 280

Max Obs per Subject 12

Optimization Information

Optimization Technique Newton-Raphson with Ridging

Parameters in Optimization 4

Lower Boundaries 3

Upper Boundaries 1

Fixed Effects Profiled

Residual Variance Profiled

Starting From Data

Iteration History

Objective Max

Iteration Restarts Evaluations Function Change Gradient

0 0 4 29327.32795 . 11583.9

1 0 2 33917.845247 -4590.517296 11739.39

Optimization routine cannot improve the function value.

Covariance Parameter Estimates

Standard

Cov Parm Subject Estimate Error

UN(1,1) ID 17.6726 .

UN(2,1) ID 1.4510 .

UN(2,2) ID 33.5211 .

SP(POW) ID 0.9000 .

Residual 1415.97 .

LOG

ERROR: NRRIDG Optimization cannot be completed.

NOTE: Optimization routine cannot improve the function value.

NOTE: PROCEDURE GLIMMIX used (Total process time):

real time 1.43 seconds

cpu time 1.37 seconds

C) CODE (Dual Quasi-Newton with random intercept)

proc glimmix data=data2 ORDER=Data ;

class id sex age xtime;

model vas = time time*time

sex age sex*age

sex*time age*time sex*age*time

sex*time*time age*time*time sex*age*time*time / ddfm=kr s cl ;

random int time time*time/ subject=id type=un v;

random _residual_ / sub=id type=sp(pow)(xtime);

run;

Model information remains same

Dimensions

G-side Cov. Parameters 6

R-side Cov. Parameters 2

Columns in X 27

Columns in Z per Subject 3

Subjects (Blocks in V) 280

Max Obs per Subject 12

Optimization Information

Optimization Technique Dual Quasi-Newton

Parameters in Optimization 7

Lower Boundaries 4

Upper Boundaries 1

Fixed Effects Profiled

Residual Variance Profiled

Starting From Data

Iteration History

Objective Max

Iteration Restarts Evaluations Function Change Gradient

0 0 4 30662.837656 . 8506.067

1 0 5 28668.856677 1993.9809793 2199.816

2 0 58 26208.493502 2460.3631744 263595.8

3 0 5 26208.485885 0.00761685 270900.4

4 0 3 26188.347277 20.13860824 308956.5

5 0 4 26158.678073 29.66920398 616998

6 0 5 26140.803057 17.87501606 654371.3

7 0 5 26124.120621 16.68243649 688251

8 0 5 26107.310639 16.80998203 722089.4

9 0 5 26090.520122 16.79051661 755560.3

10 0 5 26073.888382 16.63173957 788085.5

11 0 5 26057.541672 16.34671061 818876.5

12 0 5 26041.594314 15.94735789 846943.3

13 0 5 26026.150001 15.44431333 871095.6

14 0 5 26011.302731 14.84726934 889940.6

15 0 2 25964.073338 47.22939308 11515841

16 1 20 25867.099209 96.97412875 635971.6

17 1 3 25828.607906 38.49130370 51872.11

18 1 11 25635.995456 192.61245002 6136433

19 1 7 25632.646006 3.34945002 13163228

20 1 6 25628.79396 3.85204595 13931469

21 1 3 25605.60941 23.18454998 10127958

22 1 3 25568.842337 36.76707290 4627900

80 2 5 23953.312146 0.00000207 63380.05

81 2 4 23953.311488 0.00065810 59514.87

82 2 3 23953.311486 0.00000236 53165.53

Convergence criterion (GCONV=1E-8) satisfied.

Estimated G matrix is not positive definite.

Covariance Parameter Estimates

Standard

Cov Parm Subject Estimate Error

UN(1,1) ID 9.86E-22 .

UN(2,1) ID -3.9026 0.8652

UN(2,2) ID 2.4122 0.2092

UN(3,1) ID 0.2268 0.06188

UN(3,2) ID -0.1050 0.008000

UN(3,3) ID 7.877E-8 .

SP(POW) ID 0.09958 0.01372

Residual 140.27 4.8014

LOG

NOTE: Convergence criterion (GCONV=1E-8) satisfied.

NOTE: At least one element of the gradient is greater than 1e-3.

NOTE: Estimated G matrix is not positive definite.

NOTE: PROCEDURE GLIMMIX used (Total process time):

real time 25.57 seconds

cpu time 25.30 seconds

D) CODE (Dual Quasi-Newton without random intercept)

proc glimmix data=data2 ORDER=Data ;

class id sex age xtime;

model vas = time time*time

sex age sex*age

sex*time age*time sex*age*time

sex*time*time age*time*time sex*age*time*time / ddfm=kr s cl ;

random time time*time/ subject=id type=un v;

random _residual_ / sub=id type=sp(pow)(xtime);

run;

Model information remains same

Dimensions

G-side Cov. Parameters 3

R-side Cov. Parameters 2

Columns in X 27

Columns in Z per Subject 2

Subjects (Blocks in V) 280

Max Obs per Subject 12

Optimization Information

Optimization Technique Dual Quasi-Newton

Parameters in Optimization 4

Lower Boundaries 3

Upper Boundaries 1

Fixed Effects Profiled

Residual Variance Profiled

Starting From Data

Iteration History

Objective Max

Iteration Restarts Evaluations Function Change Gradient

0 0 4 29327.32795 . 11583.9

1 0 5 29281.411327 45.91662314 12868.42

2 0 5 29240.295101 41.11622587 14297.88

3 0 70 26042.998139 3197.2969629 5.6897E9

4 0 4 25979.595376 63.40276297 9.877E8

5 0 3 25957.25277 22.34260529 1.139E10

6 0 3 25939.317553 17.93521689 1.29E10

7 0 18 25841.391645 97.92590789 7.1394E8

8 0 11 25819.989261 21.40238439 9.1343E8

9 0 10 25767.475697 52.51356388 1.5118E9

10 0 6 25765.64274 1.83295751 2.6814E9

11 0 4 25735.316611 30.32612824 3.7442E9

12 0 2 25707.023322 28.29328991 7.3018E9

13 0 2 25665.354323 41.66899876 1.3357E9

14 0 5 25649.042009 16.31231406 5.624E9

15 0 4 25599.39337 49.64863832 1.5771E9

16 0 3 25591.215318 8.17805200 4.4921E9

17 0 4 25575.47631 15.73900814 1.106E10

18 0 2 25550.179798 25.29651182 3.693E9

19 0 5 25522.241508 27.93829080 6.2534E9

20 0 3 25496.791043 25.45046438 6.5937E9

21 0 2 25471.721861 25.06918217 2.1984E9

22 0 9 25439.671707 32.05015405 2.8937E9

130 0 2 24125.065516 22.98137596 1.1381E8

131 0 3 24112.207432 12.85808339 1.2729E8

132 0 5 24107.301953 4.90547927 55729947

133 0 3 24106.054818 1.24713534 50910488

134 0 3 24105.477952 0.57686577 9767172

135 0 3 24105.407721 0.07023051 993324.6

136 0 3 24105.40601 0.00171169 106808.8

137 0 3 24105.405963 0.00004618 89331.43

Convergence criterion (GCONV=1E-8) satisfied.

Estimated G matrix is not positive definite.

LOG

NOTE: Convergence criterion (GCONV=1E-8) satisfied.

NOTE: At least one element of the gradient is greater than 1e-3.

NOTE: Estimated G matrix is not positive definite.

NOTE: PROCEDURE GLIMMIX used (Total process time):

real time 29.81 seconds

cpu time 29.12 seconds

I hope my presentation of above facts is clear. Please let me know if i can provide any additional information.

Thanks.

Regards,

Sandhu

Re: Nonlinear Mixed Effect Model Analysis

SteveDenham — Thu, 03 Jan 2013 15:32:42 GMT

There is a reason that the default technique for optimization in GLIMMIX is the dual quasi-Newton method. It is a superior technique, and is finding an optimum solution, while the ridging method is not. You might be able to force the ridging technique by increasing the ridge parameter.

What is interesting to me is that two of the variances are essentially zero (UN(1,1) and UN(3,3)). This implies that the random effects are overspecified

Please try:

proc glimmix data=data2 ORDER=Data ;

class id sex age xtime;

model vas = time time*time

sex age sex*age

sex*time age*time sex*age*time

sex*time*time age*time*time sex*age*time*time / ddfm=kr s cl ;

random int time time*time/ subject=id type=chol v; /* Cholesky parameterization of the unstructured matrix */

random _residual_ / sub=id type=sp(pow)(xtime);

run;

If this still results in zeroes for the id and quadratic term variances, then I think some more thought needs to be given to whether the random _residual_ should be included or not. I think there is an internal contradiction between the specification of time as a continuous covariate with a random effect and a power law diminution of variation using the spatial power structure. In other words, we are overthinking the whole process. Perhaps a simpler approach will tell us all that we need.

Steve Denham

Re: Nonlinear Mixed Effect Model Analysis

Muller — Thu, 03 Jan 2013 20:23:14 GMT

Dears Steve and Sandhu,

While I do not have access to your data to make more concrete suggestions, I am not clear why you need to structure the R-side covariance matrix to the extent of what you have "type=sp(pow)(xtime)". Is the outcome variable, pain assessment using 100 mm VAS scale, measured with a lot of error that could not be accounted for using a diagonal covariance matrix as in "random _residual_/ sub=id type=VC;"? I would go with a well structured G-side covariance matrix with the default R-side (you still have a non-diagonal V matrix that accounts for both with-in and between group covariation, BTW) for both model parsimony and computational simplicity purposes. Thanks.

Mulugeta

Re: Nonlinear Mixed Effect Model Analysis

SteveDenham — Thu, 03 Jan 2013 20:35:21 GMT

Thanks for pulling us back from the brink. This is where the last two lines that I wrote was leading, but this is much more to the point, and WAAAAAY more likely to give interpretable results.

Steve Denham

Re: Nonlinear Mixed Effect Model Analysis

Sandhu — Fri, 04 Jan 2013 04:49:27 GMT

Dear Steve and Mulugeta,

I tried both the alternatives as recommended by both of you. By using Cholesky parameterization of the unstructured matrix, the only the quadratic term variances is zero, and not the id. But the use of random _residual_/ sub=id type=VC did not improve the model and it still resulted in zeroes for the id and quadratic term variances.

A) using Steve's code with Cholesky parameterization of the unstructured matrix

Covariance Parameter Estimates

Standard

Cov Parm Subject Estimate Error

CHOL(1,1) ID 3.1811 0.4606

CHOL(2,1) ID -0.1228 0.1851

CHOL(2,2) ID 0.000026 0.1909

CHOL(3,1) ID -0.00994 0.01326

CHOL(3,2) ID -2.09E-6 0.01529

CHOL(3,3) ID 4.38E-18 .

SP(POW) ID -0.1072 0.01474

Residual 115.98 3.1069

Fit Statistics

-2 Res Log Likelihood 24954.65

AIC (smaller is better) 24968.65

AICC (smaller is better) 24968.68

BIC (smaller is better) 24994.09

CAIC (smaller is better) 25001.09

HQIC (smaller is better) 24978.85

Generalized Chi-Square 388284.9

Gener. Chi-Square / DF 115.98

it converged in 49 Iteration

Log:

NOTE: Convergence criterion (GCONV=1E-8) satisfied.

NOTE: At least one element of the gradient is greater than 1e-3.

NOTE: Estimated G matrix is not positive definite.

NOTE: PROCEDURE GLIMMIX used (Total process time):

real time 25.71 seconds

cpu time 25.36 seconds

B) using Mulugeta's code with diagonal covariance matrix as in "random _residual_/ sub=id type=VC and unstructured matrix UN variance for random effect

Covariance Parameter Estimates

Standard

Cov Parm Subject Estimate Error

UN(1,1) ID 3.23E-18 .

UN(2,1) ID 0.1491 0.9892

UN(2,2) ID 1.5762 0.1915

UN(3,1) ID -0.04960 0.07431

UN(3,2) ID -0.07712 0.007403

UN(3,3) ID 0.000025 .

Residual (VC) 149.52 3.9873

Fit Statistics

-2 Res Log Likelihood 25206.89

AIC (smaller is better) 25218.89

AICC (smaller is better) 25218.92

BIC (smaller is better) 25240.70

CAIC (smaller is better) 25246.70

HQIC (smaller is better) 25227.64

it converged in 53 iterations

LOG

NOTE: Convergence criterion (GCONV=1E-8) satisfied.

NOTE: At least one element of the gradient is greater than 1e-3.

NOTE: Estimated G matrix is not positive definite.

NOTE: PROCEDURE GLIMMIX used (Total process time):

real time 17.31 seconds

cpu time 17.05 seconds

B) using Mulugeta's code with diagonal covariance matrix as in "random _residual_/ sub=id type=VC and Cholesky parameterization of the unstructured matrix variance for random effect

Covariance Parameter Estimates

Standard

Cov Parm Subject Estimate Error

CHOL(1,1) ID 0.6879 0.8774

CHOL(2,1) ID 0.1156 0.2015

CHOL(2,2) ID 4.1E-19 .

CHOL(3,1) ID -0.01065 0.01730

CHOL(3,2) ID 0.000451 0.003211

CHOL(3,3) ID 2.13E-16 .

Residual (VC) 139.60 3.6281

Fit Statistics

-2 Res Log Likelihood 26165.98

AIC (smaller is better) 26175.98

AICC (smaller is better) 26176.00

BIC (smaller is better) 26194.15

CAIC (smaller is better) 26199.15

HQIC (smaller is better) 26183.27

Generalized Chi-Square 467391.5

Gener. Chi-Square / DF 139.60

it converged in 41 iterations

Log

NOTE: Convergence criterion (GCONV=1E-8) satisfied.

NOTE: At least one element of the gradient is greater than 1e-3.

NOTE: Estimated G matrix is not positive definite.

NOTE: PROCEDURE GLIMMIX used (Total process time):

real time 3.25 seconds

cpu time 3.19 seconds

Thanks

Regards

Sandhu

Re: Nonlinear Mixed Effect Model Analysis

SteveDenham — Fri, 04 Jan 2013 12:38:55 GMT

It is no great crime if a variance component is zero (G matrix is not positive definite). It simply means that for the data available, other terms explain all of the variation. So long as the proper adjustment to degrees of freedom is applied, tests and confidence intervals should be acceptable.

Given all that has gone on here, it appears that the Cholesky parameterization leads to the model with the smallest information criteria, so I would just go with that. It implies that after the individual variability, linear in time variability, their covariance, and their covariance with quadratic in time are accounted for, you have no remaining variability due strictly to the quadratic in time component. So long as you feel that the covariances are important, then this is where you end up. With a larger or different data set, you might get something estimated for the quadratic component.

If you believe, a priori, that the covariances are not important, then you should steer away from the unstructured covariance types, and use the diagonal structure that Muller proposed.

Steve Denham

Re: Nonlinear Mixed Effect Model Analysis

Sandhu — Fri, 04 Jan 2013 12:52:39 GMT

Thanks a lot, Steve. I will do that.

Regrads

Sandhu

Re: Nonlinear Mixed Effect Model Analysis

Muller — Fri, 04 Jan 2013 14:02:29 GMT

I agree with Steve. Also, you need to separate modeling the mean part (what terms of time to include... linear, quadratic, etc) and modeling the covariance part. The book by Fitzmaurice et al has a good discussion of steps to follow when modeling the mean and covariance components of GLMM. I usually like to identify the 'best' mean model under a more unstructured covariance model and then work towards a more parsimonious covariance model using my 'best' mean model. Hope this helps.

Mulugeta

Re: Nonlinear Mixed Effect Model Analysis

Sandhu — Sat, 05 Jan 2013 13:04:17 GMT

Dear Mulugeta,

Are you referring to "Longitudnal data analysis 2009" book by Fitzmaurice et al? i tried very hard to get this book but could not. I really feel this is very good book as per little bit information i could get from this book through internet. I will follow your's and Steve's advise for analysis.

Thanks a lot for all the help.

Sandhu

Re: Nonlinear Mixed Effect Model Analysis

jrbrauer — Tue, 15 Jan 2013 16:02:33 GMT

Hi all,

I'm sorry for such a late addition to the discussion here - just ignore if you have solved this problem already. In all the examples provided, you appear to start with the assumption that a baseline quadratic growth model offers the best fit for the data. You might consider revisiting this assumption. The only post that provides any fixed effects coefficients (#6) suggests that a linear model might fit as well (or better). At the very least, the fixed portion of the model appears seriously overparameterized; further, inclusion of squared terms for time (implying a common curvature in VAS scores across individuals) and interactions with the squared term (implying systematic differences in temporal curvature across age and sex groups) in this portion of the model might be leading you to overparameterize the stochastic portion of the model as well (as Steve noted earlier). Have you tried hierarchically building the model from a simpler platform, perhaps like the following:

*model baseline linear growth model;

proc glimmix data=data2 ORDER=Data ;

class id sex age ;

model vas = int time;

random int time / subject=id type=un;

run;

*model baseline quadratic growth model and compare overall fit and statistical significance of squared term;

proc glimmix data=data2 ORDER=Data ;

class id sex age ;

model vas = int time time*time

random int time / subject=id type=un;

run;

*model baseline quadratic growth model with additional error term, compare overall fit and statistical significance of error term;

proc glimmix data=data2 ORDER=Data ;

class id sex age ;

model vas = int time time*time

random int time time*time / subject=id type=un;

run;

*from what I've seen so far, at this point, I'm guessing you will need to remove time*time from the stochastic portion of the model, and perhaps remove its fixed coefficient as well;

*for instance, you may be moving to something like the following...;

*model average effect of sex, age, and their interaction on the intercept (mean at time 0) of VAS trajectories

proc glimmix data=data2 ORDER=Data ;

class id sex age ;

model vas = int time time*time

sex age sex*age

random int time / subject=id type=un;

run;

*then, continue to build on the model to determine whether additional interactions are warranted, for instance... ;

*model average effects of sex, age, and their interaction on the intercept (mean at time 0) and linear slope (instantaneous rate of change) of VAS trajectories;

proc glimmix data=data2 ORDER=Data ;

class id sex age ;

model vas = int time time*time

sex age sex*age

sex*time age*time sex*age*time

random int time / subject=id type=un;

run;

*again, based on predictions from post #6, I'm guessing additional interactions (and perhaps even these) are unwarranted (i.e., age*time*time, etc.);

After building the fixed portion of the model, you might then attempt to improve model fit by further specifying the stochastic portion (I think this approach is akin to what Muller was recently suggesting). For instance, you might then consider specifying an autoregressive covariance structure to account for stronger correlations between proximate observations in a repeated measures design. I agree with Muller that an R-side correction is probably unnecessary, but at this point, you could compare...

I hope this helps and that it is not too late.

Best,

Jon

Message was edited by: Jon Brauer

Re: Nonlinear Mixed Effect Model Analysis

Sandhu — Wed, 16 Jan 2013 14:25:49 GMT

Hello Jon,

Thanks for such a nice explanation. I will consider these options also. Its great to have such a helpful discussion with you, Steve and Muller.

Thanks

Sandhu