Simulations in SAS

Gompertz - 2

%let nsim=100;
data power;
length TT $30;
array a[4] _temporary_ (314.31 		347 		314.31 		308.4);
array k[4] _temporary_ (0.06158 	0.00975 	0.06158 	0.0296);
array t[4] _temporary_ (-15.442 	-205.4 		-18.9413 	-63.43); /* You can already see that the estimates of T are not precies and helpfull --> carefull attention is needed for these */
do sim=1 to ≁
	do treatment=1 to 4;
	if treatment=1 then TT="A";
	if treatment=2 then TT="B";
	if treatment=3 then TT="C";
	if treatment=4 then TT="D";
		do dose=0,20,40,80,160;
			do pen=1 to 6;
				if TT="A" 	then Y=a[1]*exp(-exp(-k[1]*(dose-t[1]))) + rand("Normal", 0, sqrt(480.779));
				if TT="B" 	then Y=a[2]*exp(-exp(-k[2]*(dose-t[2]))) + rand("Normal", 0, sqrt(480.779));
				if TT="C" 	then Y=a[3]*exp(-exp(-k[3]*(dose-t[3]))) + rand("Normal", 0, sqrt(480.779));
				if TT="D" 	then Y=a[4]*exp(-exp(-k[4]*(dose-t[4]))) + rand("Normal", 0, sqrt(480.779));		/* What is the true variance!!! */
				output;
			end;
		end;
	end;
end;
run;
proc sort data=power;by sim TT dose;run;
/* Check Simulation via Graphs */
ods graphics / width=10in;
proc sgpanel data=power;
panelby sim / columns=5;
where sim<11;
vbox Y / category=dose group=TT;
run;
ods graphics / reset;
/* 
 * Check Simulation via NLIN
 */
proc nlin data=power plots=all maxiter=20 hougaard nlinmeasures bias;
where sim=1 and treatment=1;
parms a=314 k=0.06 T=-15;
model Y=a*exp(-exp(-k*(Dose-T)));
output out=nlinpred p=yhat r=resid student=stresid lcl=lpl lclm=lcl ucl=upl uclm=ucl parms=a k T;
run;
proc sgplot data=nlinpred;
scatter y=Y   x=dose;
series  y=yhat   	x=dose / name="fit" lineattrs=(color=red pattern=2) LegendLabel="Predicted";
series  y=lpl 		x=dose / lineattrs=GraphConfidence name="lpl" legendlabel="95% Prediction limits";
series  y=upl 		x=dose / lineattrs=GraphConfidence;
band  	x=dose lower=lcl upper=ucl / fillattrs=(color=red) name="lcl" legendlabel="95% Confidence limits" transparency=0.75;
refline k / axis=x lineattrs=(color=orange pattern=2) label="Inflection Point";
refline a / axis=y lineattrs=(color=green pattern=2) label="Asymptot";
keylegend "fit" "lpl" "lcl";
yaxis min=0 max=350 label="Y  ";
xaxis label="Day  ";
Title "Gompertz to describe Y   ";
run;
proc nlin data=power plots=all maxiter=20 hougaard nlinmeasures bias;
where sim=1 and treatment=3;
parms a=314 k=0.06 T=-15;
model Y=a*exp(-exp(-k*(Dose-T)));
output out=nlinpred p=yhat r=resid student=stresid lcl=lpl lclm=lcl ucl=upl uclm=ucl parms=a k T;
run;
proc sgplot data=nlinpred;
scatter y=Y   		x=dose;
series  y=yhat   	x=dose / name="fit" lineattrs=(color=red pattern=2) LegendLabel="Predicted";
series  y=lpl 		x=dose / lineattrs=GraphConfidence name="lpl" legendlabel="95% Prediction limits";
series  y=upl 		x=dose / lineattrs=GraphConfidence;
band  	x=dose lower=lcl upper=ucl / fillattrs=(color=red) name="lcl" legendlabel="95% Confidence limits" transparency=0.75;
refline k / axis=x lineattrs=(color=orange pattern=2) label="Inflection Point";
refline a / axis=y lineattrs=(color=green pattern=2) label="Asymptot";
keylegend "fit" "lpl" "lcl";
yaxis min=0 max=350 label="Y  ";
xaxis label="Day  ";
Title "Gompertz to describe Y   ";
run;