Re: PROC MIXED

GiacomoF · Posted 01-20-2016 06:58 AM

Suppose to have a Dataset with the following variables Simulation (1 to 1000),Study (1 to 5), Treated (Dummy variable 0 not treated, 1 treated), Outcome (1 to 3), Outcome1 (Dummy variable 0 is not measured, 1 is measured), Outcome 2 (the same of Outcome1), Outcome3 (the same of Outcome1 and Outcome2), ID: identification of patient, Y: response variable.
I has used proc mixed written as follow:

proc mixed cl method=reml data=....;
by simulation study;
class Outcome Id Treated;
model Y=outcome1 outcome2 outcome3 treated*outcome1 treated*outcome2 treated*outcome3 /noint s cl covb;
repeated Outcome/ type=un subject=id;
run;

to obtain the fixed effects estimates that are the effects for each outcome across the groups of treated and not treated patients.

If I wish to obtain the random effects for the estimates of the treatment across the groups for each simulation and study how should I change the proc mixed I wrote before?

proc mixed cl method=reml data=....;
by simulation study;
class Outcome Id Treated;
model Y= /;

random outcome1 outcome2 outcome3 treated*outcome1 treated*outcome2 treated*outcome3/s cl;
run;

If i run this I obtain to much DF

proc mixed cl method=reml data=....;
by simulation study;
class Outcome Id Treated;
model Y= /;

random outcome1 outcome2 outcome3 treated*outcome1 treated*outcome2 treated*outcome3/s cl type=un subject=id;
run;

If i run this The hessian matrix is not positive definite therefore I not obtain any estimates.

Do you understand the problem.....I have also a computational problems because if I change the code considering the command RANDOM in the proc mixed SAS seems to take a lot of minuts.

Rick_SAS · Posted 01-20-2016 08:13 AM

Often "hessian not positive definite" occurs when the model is not appropriate for the data or the sample size is too small.

What is your sample size? Do you get the error for all simulations, or just some?

I am confused by your models. When you simulate the data, what is the model for the simulation? Are the TREATED and OUTSOME: variables simulated as fixed or random effects? I recommend that you post your simulation code so we can see what you are simulating.

GiacomoF · Posted 01-20-2016 08:24 AM

To simulate I use the following proc iml. The aim is to simulate random effects......

proc iml;

N1t = 30; N2t = 40; N3t = 50; N4t = 60; N5t = 70;

N1p = 30; N2p = 40; N3p = 50; N4p = 60; N5p = 70;

NumSamples=1000;

MeanT = {4, 5, 15};

MeanP = {9, 8, 21};

CovT = {74 0 0, 0 80 0, 0 0 117};

CovP = {49 0 0, 0 64 0, 0 0 81};

call randseed(100);

ST1 = RandNormal(N1t*NumSamples, MeanT, CovT);

ST2 = RandNormal(N2t*NumSamples, MeanT, CovT);

ST3 = RandNormal(N3t*NumSamples, MeanT, CovT);

ST4 = RandNormal(N4t*NumSamples, MeanT, CovT);

ST5 = RandNormal(N5t*NumSamples, MeanT, CovT);

SP1 = RandNormal(N1p*NumSamples, MeanP, CovP);

SP2 = RandNormal(N2p*NumSamples, MeanP, CovP);

SP3 = RandNormal(N3p*NumSamples, MeanP, CovP);

SP4 = RandNormal(N4p*NumSamples, MeanP, CovP);

SP5 = RandNormal(N5p*NumSamples, MeanP, CovP);

S1t = colvec(repeat(T(1:NumSamples), 1, N1t));

S2t = colvec(repeat(T(1:NumSamples), 1, N2t));

S3t = colvec(repeat(T(1:NumSamples), 1, N3t));

S4t = colvec(repeat(T(1:NumSamples), 1, N4t));

S5t = colvec(repeat(T(1:NumSamples), 1, N5t));

S1p = colvec(repeat(T(1:NumSamples), 1, N1p));

S2p = colvec(repeat(T(1:NumSamples), 1, N2p));

S3p = colvec(repeat(T(1:NumSamples), 1, N3p));

S4p = colvec(repeat(T(1:NumSamples), 1, N4p));

S5p = colvec(repeat(T(1:NumSamples), 1, N5p));

Z1t = S1t || ST1; Z2t = S2t || ST2; Z3t = S3t || ST3; Z4t = S4t || ST4; Z5t = S5t || ST5;

Z1p = S1p || SP1; Z2p = S2p || SP2; Z3p = S3p || SP3; Z4p = S4p || SP4; Z5p = S5p || SP5;

create work.Treated1 from Z1t[c={"Simulation" "x1" "x2" "x3"}]; append from Z1t;

close work.Treated1;

create work.Treated2 from Z2t[c={"Simulation" "x1" "x2" "x3"}]; append from Z2t;

close work.Treated2;

create work.Treated3 from Z3t[c={"Simulation" "x1" "x2" "x3"}]; append from Z3t;

close work.Treated3;

create work.Treated4 from Z4t[c={"Simulation" "x1" "x2" "x3"}]; append from Z4t;

close work.Treated4;

create work.Treated5 from Z5t[c={"Simulation" "x1" "x2" "x3"}]; append from Z5t;

close work.Treated5;

create work.Placebo1 from Z1p[c={"Simulation" "x1" "x2" "x3"}]; append from Z1p;

close work.Placebo1;

create work.Placebo2 from Z2p[c={"Simulation" "x1" "x2" "x3"}]; append from Z2p;

close work.Placebo2;

create work.Placebo3 from Z3p[c={"Simulation" "x1" "x2" "x3"}]; append from Z3p;

close work.Placebo3;

create work.Placebo4 from Z4p[c={"Simulation" "x1" "x2" "x3"}]; append from Z4p;

close work.Placebo4;

create work.Placebo5 from Z5p[c={"Simulation" "x1" "x2" "x3"}]; append from Z5p;

close work.Placebo5;

quit;

Rick_SAS · Posted 01-20-2016 09:09 AM

Your IML program doesn't seem related to your original question. The code you've supplied generates variables x1-x3 as uncorrelated normal variables. Your question involves binary variables named TREATED and OUTCOME variables. Your program also does not specify the relationship between the reponse variable Y and the explanatory variables.

GiacomoF · Posted 01-20-2016 09:27 AM

x1 is the response variable for outcome 1

x2 is the response variable for outcome 2

x3 is the response variable for outcome 3

You are right when you say that are uncorrelated normal variables as in the Covariance Matrix I assumes 0 correlation.

You are also right that this IML code does not consider the variable Treated, Outcome and id.

/*To introduce Treated (0 is not treated; 1 is treated)*/

data work.Treated1; set work.Treated1; Study=1; Treated=0; run;
data work.Treated2; set work.Treated2; Study=2; Treated=0; run;
data work.Treated3; set work.Treated3; Study=3; Treated=0; run;
data work.Treated4; set work.Treated4; Study=4; Treated=0; run;
data work.Treated5; set work.Treated5; Study=5; Treated=0; run;
data work.Placebo1; set work.Placebo1; Study=1; Treated=1; run;
data work.Placebo2; set work.Placebo2; Study=2; Treated=1; run;
data work.Placebo3; set work.Placebo3; Study=3; Treated=1; run;
data work.Placebo4; set work.Placebo4; Study=4; Treated=1; run;
data work.Placebo5; set work.Placebo5; Study=5; Treated=1; run;

/*To introduce the Id for each patients*/

data work.idT1; retain Seed1 12985062; do simulation=1 to 1000; do id=1 to 30;
X1=12.1*rannor(Seed1)+11.7; output; end; end; run;
data work.idT2; retain Seed1 12985062; do simulation=1 to 1000; do id=1 to 40;
X1=12.1*rannor(Seed1)+11.7; output; end; end; run;
data work.idT3; retain Seed1 12985062; do simulation=1 to 1000; do id=1 to 50;
X1=12.1*rannor(Seed1)+11.7; output; end; end; run;
data work.idT4; retain Seed1 12985062; do simulation=1 to 1000; do id=1 to 60;
X1=12.1*rannor(Seed1)+11.7; output; end; end; run;
data work.idT5; retain Seed1 12985062; do simulation=1 to 1000; do id=1 to 70;
X1=12.1*rannor(Seed1)+11.7; output; end; end; run;
data work.idP1; retain Seed1 12985062; do simulation=1 to 1000; do id=31 to 60;
X1=12.1*rannor(Seed1)+11.7; output; end; end; run;
data work.idP2; retain Seed1 12985062; do simulation=1 to 1000; do id=41 to 80;
X1=12.1*rannor(Seed1)+11.7; output; end; end; run;
data work.idP3; retain Seed1 12985062; do simulation=1 to 1000; do id=51 to 100;
X1=12.1*rannor(Seed1)+11.7; output; end; end; run;
data work.idP4; retain Seed1 12985062; do simulation=1 to 1000; do id=61 to 120;
X1=12.1*rannor(Seed1)+11.7; output; end; end; run;
data work.idP5; retain Seed1 12985062; do simulation=1 to 1000; do id=71 to 140;
X1=12.1*rannor(Seed1)+11.7; output; end; end; run;

/*I keep simulation Id*/
data work.idT1; set work.idT1; keep Simulation Id; run;
data work.idT2; set work.idT2; keep Simulation Id; run;
data work.idT3; set work.idT3; keep Simulation Id; run;
data work.idT4; set work.idT4; keep Simulation Id; run;
data work.idT5; set work.idT5; keep Simulation Id; run;
data work.idP1; set work.idP1; keep Simulation Id; run;
data work.idP2; set work.idP2; keep Simulation Id; run;
data work.idP3; set work.idP3; keep Simulation Id; run;
data work.idP4; set work.idP4; keep Simulation Id; run;
data work.idP5; set work.idP5; keep Simulation Id; run;

/*I merge the datasets for each studies to the column Id*/
data work.T1; merge work.Treated1 work.idT1; run;
data work.T2; merge work.Treated2 work.idT2; run;
data work.T3; merge work.Treated3 work.idT3; run;
data work.T4; merge work.Treated4 work.idT4; run;
data work.T5; merge work.Treated5 work.idT5; run;
data work.P1; merge work.Placebo1 work.idP1; run;
data work.P2; merge work.Placebo2 work.idP2; run;
data work.P3; merge work.Placebo3 work.idP3; run;
data work.P4; merge work.Placebo4 work.idP4; run;
data work.P5; merge work.Placebo5 work.idP5; run;

/*I create one dataset for each Study*/

proc datasets; append base=work.Study1 data=work.T1; run;
proc datasets; append base=work.Study1 data=work.P1; run;
proc datasets; append base=work.Study2 data=work.T2; run;
proc datasets; append base=work.Study2 data=work.P2; run;
proc datasets; append base=work.Study3 data=work.T3; run;
proc datasets; append base=work.Study3 data=work.P3; run;
proc datasets; append base=work.Study4 data=work.T4; run;
proc datasets; append base=work.Study4 data=work.P4; run;
proc datasets; append base=work.Study5 data=work.T5; run;
proc datasets; append base=work.Study5 data=work.P5; run;

/*I create one dataset called Datasim1*/
proc datasets; append base=work.Datasim1 data=work.Study1; run;
proc datasets; append base=work.Datasim1 data=work.Study2; run;
proc datasets; append base=work.Datasim1 data=work.Study3; run;
proc datasets; append base=work.Datasim1 data=work.Study4; run;
proc datasets; append base=work.Datasim1 data=work.Study5; run;

/*I create the dummy variable Outcome 1*/

data work.Outcome1 (keep=Simulation Study Id outcome1 outcome2 outcome3 Treated Outcome Y);
set work.Datasim1;
outcome1=1; outcome2=0; outcome3=0; Outcome=1; Y=x1;
run;

/*I create the dummy variable Outcome 2*/
data work.Outcome2 (keep=Simulation Study Id outcome1 outcome2 outcome3 Treated Outcome Y);
set work.Datasim1;
outcome1=0; outcome2=1; outcome3=0; Outcome=2; Y=x2;
run;

/*I create the dummy variable Outcome 3*/
data work.Outcome3 (keep=Simulation Study Id outcome1 outcome2 outcome3 Treated Outcome Y);
set work.Datasim1;
outcome1=0; outcome2=0; outcome3=1; Outcome=3; Y=x3;
run;

/*I create the Dataset that contains Dummy variables and Y response variable*/

proc datasets; append base=work.Datasim2 data=work.Outcome1; run;
proc datasets; append base=work.Datasim2 data=work.Outcome2; run;
proc datasets; append base=work.Datasim2 data=work.Outcome3; run;

proc sort data=work.Datasim2;
by Simulation Study Outcome;
run;

Rick_SAS · Posted 01-20-2016 10:08 AM

I encourage you to take a step back and to think about the data that you are obtaining for each simulation and study. For example, restrict your attention to the data set where simulation=1 and study=1. Given that ONE data set, figure our what statistical questions you are trying to answer.

For study=1, it looks like you are trying to simulate two groups of 30 patients. One group is treated. The other is untreated. All patients come into a clinic three times and have an outcome measured. You can visualize this situation by using a panel of spaghetti plots:

data A;
set Datasim2;
where simulation=1 AND study=1;
run;

title "Simulation=1; Study=1";
proc sgpanel data=A;
panelby treated;
series x=outcome y=Y / group=id;
run;

Is this correct? Does this represent the real-life scenario that you are trying to simulate?

After you've shown an example and explained what you are trying to measure and compare (mean slope?), the mixed-model experts in thei community can help you determine what statistical questions to ask and what PROC MIXED code you need to contruct. Good luck!

lvm · Posted 01-20-2016 12:30 PM

Rick makes good points. I don't see anything in your simulation to give you the mixed model you are trying to fit. You are trying to fit a very complex variance-covariance model. I can see why you would have problems with overparameterization.

GiacomoF · Posted 01-21-2016 08:34 AM

To answer to IVM when I use the command RANDOM in the proc mixed I have problem of DF. To solve this problem I think the solution is to say to the proc mixed the correct number of observation to use for each simulation and study......because for example for study 5 where I have 140 subjects the DF calculated form the proc mixed model are (140*3)-6=414.

proc mixed cl method=reml data=work.Datasim2;
by Simulation Study;
class Outcome id Treated;
ods listing close;
model Y = / noint;
random Outcome1 Outcome2 Outcome3 Treated*Outcome1 Treated*Outcome2 Treated*Outcome3/ s cl;
run;

SteveDenham · Posted 01-22-2016 10:04 AM

Random effects don't use up degrees of freedom in SAS mixed models the same way that they do in ordinary ANOVA (depends on the covariance structure chosen). If the Kenward-Rogers method is giving you unusual results, then consider using ddfm=bw (between-within) and some sort of empirical shrinkage.

Steve Denham

GiacomoF · Posted 01-22-2016 10:08 AM

Hello Steve DDFM in the random statment?

SteveDenham · Posted 01-22-2016 10:19 AM

In the MODEL statement. The RANDOM statement will not take a degrees of freedom modifier, as the estimates are asymptotes to the true variance components which are N(0,<value>) distributed. No df's involved.

Steve Denham

GiacomoF · Posted 01-21-2016 08:23 AM

To answer to Rick, yes I tried to look at the graphics as you suggested and I obtain the same scenario of the graphics you posted yesterday.