Zhang, Z., & Wang, L. (2009). Power analysis for growth curve models using SAS. Behavior Research Methods, 41(4), 1083-1094. Request a copy
/*Suppress the output and the log */
options nosource nonotes nosource2 nomprint;
/*CHANGE THE PARAMETERS HERE*/
*model parameters;
%LET MuL=10; *mean peak level;
%LET MuS=5; *mean change between the initial status and the peak level;
%LET Sigma_e=1; *residual standard deviation;
%LET Sigma_L=2; *level standard deviation;
%LET Sigma_S=.5; *change standard deviation;
%LET rho=0.5; *correlation between Li and Si;
%LET miss=0; *missing data rate, 0: no missing data;
%LET p=1; *rate of growth/decline;
*power parameters;
%LET R=10000; *number of simulation replications;
%LET T=6; *number of measurement occasions;
%LET start=100; *the minimum sample size to consider;
%LET end=100; *the maximum sample size to consider;
%LET step=50; *the step between two sample sizes;
%LET df=2; *the difference in the number of parameters;
%LET seed=1; *random number generator seed;
/*DO NOT CHANGE CODES BELOW UNLESS YOU KNOW WHAT YOU ARE DOING*/
/*Calculate the chi-square difference between two nested growth curve models*/
%MACRO LL(N,T,seed);
DATA Sim_ExpGM;
* set statistical parameters;
N = &N; seed = &seed;
* setup arrays for repeated measures;
ARRAY y_score{&T} y1-y&T;
ARRAY M{&T} m1-m&T;
m1=1;
* generate raw data with considering the missing data rate;
DO _N_ = 1 TO N;
e_L=RANNOR(seed);
e_S=&rho*e_L+SQRT(1-&rho**2)*RANNOR(seed);
L_score=&MuL+&Sigma_L*e_L;
S_score=&MuS+&Sigma_S*e_S;
* include indicator variables to generate missing data ;
DO t = 1 TO &T;
y_score{t} = L_score - S_score*exp(-(t-1)*&p) + &Sigma_e*RANNOR(seed);
END;
DO t=2 TO &T;
m{t}=m{t-1};
IF m{t-1}=1 AND RANUNI(seed) > (1-&miss * (t-1))/(1-&miss * (t-2)) THEN m{t} = m{t-1}*0;
IF m{t}=0 THEN y_score{t}=.;
END;
KEEP y1-y&T;
OUTPUT;
END;
RUN;
DATA ExpGM;
SET Sim_ExpGM;
%DO t = 1 %TO &T;
id = _N_; time=&t-1; y=y&t; OUTPUT;
%END;
KEEP id time y;
RUN;
/*Fit two nested models to the data*/
ODS OUTPUT FitStatistics(persist=proc)=fit;
*Model 1: the true model - exponential growth curve model;
PROC NLMIXED DATA = ExpGM;
traject = level-slope*exp(-p*(time));
MODEL y ~ NORMAL(traject, v_e);
RANDOM level slope ~ NORMAL([m_l, m_s], [v_l, c_ls, v_s])
SUBJECT = id;
PARMS m_l = 10 m_s=5 v_l = 4 c_ls = 0.5 v_s = .25 v_e = 1 p=1;
RUN;
*Model 2: the null model - no variation in Si;
PROC NLMIXED DATA = ExpGM;
traject = level-m_s*exp(-p*(time));
MODEL y ~ NORMAL(traject, v_e);
RANDOM level ~ NORMAL(m_l, v_l)
SUBJECT = id;
PARMS m_l = 10 m_s=5 v_l = 4 v_e = 1 p=1;
RUN;
ODS OUTPUT CLOSE;
%MEND LL;
/*The second Macro: POWER*/
/*This Macro calls the first Macro LL for each replication*/
* Calculate power based on R replications;
%MACRO POWER(R,N,T,seed,df);
DATA tempfit;
DO _N_=1 TO 8;
tempfit=_N_;
OUTPUT;
END;
RUN;
%LL(&N,&T,&seed);
DATA fit;
MERGE fit tempfit;
RUN;
DATA allfit;
SET fit;
RUN;
%DO I = 2 %TO &R;
PROC DATASETS LIBRARY=WORK; DELETE fit; RUN; QUIT;
%LL(&N,&T,%eval(&seed+&I*1389));
DATA fit;
MERGE fit tempfit;
RUN;
DATA allfit;
SET allfit fit;
RUN;
DM 'CLEAR LOG';
%END;
DATA allfit;
SET allfit;
IF MOD(_N_,4) ~= 1 THEN DELETE;
KEEP Value;
RUN;
DATA allfit;
SET allfit;
id =INT((_N_-.1)/2)+1;
modelnum = MOD(_N_+1, 2);
RUN;
PROC TRANSPOSE DATA=allfit OUT=allfit prefix=model;
BY id;
ID modelnum;
VAR Value;
RUN;
DATA allfit;
SET allfit;
ss = &N;
diff = model1 - model0;
ind = 1;
IF diff=. THEN DELETE;
IF diff<0 THEN DELETE;
IF diff < CINV(.95, &df) THEN ind = 0;
DROP id _NAME_ model0 model1;
RUN;
PROC MEANS DATA = allfit;
VAR ss ind;
OUTPUT OUT=power mean(ss ind)=ss power;
RUN;
%MEND POWER;
*%POWER(&R, 100, &T, &seed, &df);
/*The third Macro: POWERCURVE*/
/* This Macro calls the second Macro for each sample size*/
%MACRO powercurve(R, seed, st, end, step, T,df);
%POWER(&R, &st, &T, %eval(&seed+&st), &df);
DATA allpower;
SET power;
RUN;
%LET st = %eval(&st + &step);
%DO %WHILE (&st <= &end);
%POWER(&R, &st, &T, %eval(&seed+&st), &df);
%LET st = %eval(&st + &step);
DATA allpower;
SET allpower power;
RUN;
DM 'CLEAR LOG';
%END;
* Save the results for possible future use;
DATA allpower;
SET allpower;
FILE "power-exp-new.txt";
PUT _FREQ_ ss power;
RUN;
* Plot the power curve;
ODS PDF FILE='power.pdf' NOTOC;
PROC GPLOT DATA = allpower;
SYMBOL I=JOIN;
PLOT power*ss;
RUN;
QUIT;
ODS PDF CLOSE;
%MEND powercurve;
ODS RESULTS OFF;
ODS LISTING CLOSE;
%powercurve(&R,&seed,&start,&end,&step,&T,&df);
ODS RESULTS ON;
ODS LISTING;