Solved: Compute standard errors without using PROC IML

tka726 · Posted 11-12-2020 12:49 PM

Hello!

I need to use a macro which requires use of IML - which I don't have. Is there another way to program this to get the result?

For reference, the full macro is here:

https://www.sas.com/content/dam/SAS/support/en/sas-global-forum-proceedings/2018/2424-2018.pdf

/*Fit Hamlett (2004) mixed model*/
proc mixed data=&data._long method=ml covtest asycov asycorr;
class &ID. vtype &rep.;
model response = vtype / solution ddfm=kr;
random vtype / type=un subject=&ID. v vcorr;
repeated vtype / type=un subject=&rep.(&ID.);
ods output VCorr=VCorr ConvergenceStatus=CS asycorr=asycorr
asycov=asycov CovParms=CovParms;
run;

/*Use delta method to find SE(rho), calculate CI/p-value based on normality*/
proc iml;
use Covparms;
read all var {Estimate} into Cov_estimate;
close Covparms;
use Asycov;
read all var {CovP1 CovP2 CovP3 CovP4 CovP5 CovP6} into asycov;
close Asycov;
a = Cov_estimate[1];
b = Cov_estimate[2];
c = Cov_estimate[3];
g = Cov_estimate[4];
h = Cov_estimate[5];
i = Cov_estimate[6];
rhohat = (b+h)/sqrt((a+g)*(c+i));
df_da = -0.5*((b+h)*(c+i))/sqrt(((a+g)*(c+i))**3);
df_db = 1/sqrt((a+g)*(c+i));
df_dc = -0.5*(b+h)*(a+g)/sqrt(((a+g)*(c+i))**3);
df_dg = -0.5*(b+h)*(c+i)/sqrt(((a+g)*(c+i))**3);
df_dh = 1/sqrt((a+g)*(c+i));
df_di = -0.5*(b+h)*(a+g)/sqrt(((a+g)*(c+i))**3);
create partialderiv var {df_da df_db df_dc df_dg df_dh df_di};
append;
close partialderiv;
use partialderiv;
read all var {df_da df_db df_dc df_dg df_dh df_di} into partialderiv;
close partialderiv;
Sigma = asycov;
rho_var = partialderiv*Sigma*t(partialderiv);
rho_SE = sqrt(rho_var);
l95=rhohat-1.96*rho_SE;
u95=rhohat+1.96*rho_SE;
p = (1-cdf("normal",abs(rhohat),0,rho_SE))*2;
NormalApproxOutput = j(1,5,.);
NormalApproxOutput[,1]=rhohat;
NormalApproxOutput[,2]=rho_SE;
NormalApproxOutput[,3]=l95;
NormalApproxOutput[,4]=u95;
NormalApproxOutput[,5]=p;
Outtitle={'Estimate' 'Std Error' '95% CI Lower Bound' '95% CI Upper
Bound' 'Pvalue'};
Varnames=t("&var1. and &var2.");
PRINT NormalApproxOutput[colname=Outtitle rowname=Varnames];
quit;

Ksharp · Posted 11-15-2020 06:20 AM

The IML is similar with data step.
The following code could translate into date step like:
proc iml;
use Covparms;
read all var {Estimate} into Cov_estimate;
close Covparms;
use Asycov;
read all var {CovP1 CovP2 CovP3 CovP4 CovP5 CovP6} into asycov;
close Asycov;
a = Cov_estimate[1];
b = Cov_estimate[2];
c = Cov_estimate[3];
g = Cov_estimate[4];
h = Cov_estimate[5];
i = Cov_estimate[6];
rhohat = (b+h)/sqrt((a+g)*(c+i));
df_da = -0.5*((b+h)*(c+i))/sqrt(((a+g)*(c+i))**3);
df_db = 1/sqrt((a+g)*(c+i));
df_dc = -0.5*(b+h)*(a+g)/sqrt(((a+g)*(c+i))**3);
df_dg = -0.5*(b+h)*(c+i)/sqrt(((a+g)*(c+i))**3);
df_dh = 1/sqrt((a+g)*(c+i));
df_di = -0.5*(b+h)*(a+g)/sqrt(((a+g)*(c+i))**3);

---------->>

data partialderiv ;
set Covparms end=last;
retain a b d g h i;
if _n_=1 then a=Estimate;
if _n_=2 then b=Estimate;
if _n_=3 then c=Estimate;
if _n_=4 then g=Estimate;
if _n_=5 then h=Estimate;
if _n_=6 then i=Estimate;

if last then do;
rhohat = (b+h)/sqrt((a+g)*(c+i));
df_da = -0.5*((b+h)*(c+i))/sqrt(((a+g)*(c+i))**3);
df_db = 1/sqrt((a+g)*(c+i));
df_dc = -0.5*(b+h)*(a+g)/sqrt(((a+g)*(c+i))**3);
df_dg = -0.5*(b+h)*(c+i)/sqrt(((a+g)*(c+i))**3);
df_dh = 1/sqrt((a+g)*(c+i));
df_di = -0.5*(b+h)*(a+g)/sqrt(((a+g)*(c+i))**3);

output;
end;
keep df_da df_db df_dc df_dg df_dh df_di;
run;

View solution in original post

Ksharp · Posted 11-13-2020 07:14 AM

IML questions for IML Forum , @Rick_SAS is there .

tka726 · Posted 11-14-2020 08:57 AM

Hello!

I need to use a macro which requires use of IML - which I don't have. Is there another way to program this to get the result?

For reference, the full macro is here:

https://www.sas.com/content/dam/SAS/support/en/sas-global-forum-proceedings/2018/2424-2018.pdf

/*Fit Hamlett (2004) mixed model*/
proc mixed data=&data._long method=ml covtest asycov asycorr;
class &ID. vtype &rep.;
model response = vtype / solution ddfm=kr;
random vtype / type=un subject=&ID. v vcorr;
repeated vtype / type=un subject=&rep.(&ID.);
ods output VCorr=VCorr ConvergenceStatus=CS asycorr=asycorr
asycov=asycov CovParms=CovParms;
run;

/*Use delta method to find SE(rho), calculate CI/p-value based on normality*/
proc iml;
use Covparms;
read all var {Estimate} into Cov_estimate;
close Covparms;
use Asycov;
read all var {CovP1 CovP2 CovP3 CovP4 CovP5 CovP6} into asycov;
close Asycov;
a = Cov_estimate[1];
b = Cov_estimate[2];
c = Cov_estimate[3];
g = Cov_estimate[4];
h = Cov_estimate[5];
i = Cov_estimate[6];
rhohat = (b+h)/sqrt((a+g)*(c+i));
df_da = -0.5*((b+h)*(c+i))/sqrt(((a+g)*(c+i))**3);
df_db = 1/sqrt((a+g)*(c+i));
df_dc = -0.5*(b+h)*(a+g)/sqrt(((a+g)*(c+i))**3);
df_dg = -0.5*(b+h)*(c+i)/sqrt(((a+g)*(c+i))**3);
df_dh = 1/sqrt((a+g)*(c+i));
df_di = -0.5*(b+h)*(a+g)/sqrt(((a+g)*(c+i))**3);
create partialderiv var {df_da df_db df_dc df_dg df_dh df_di};
append;
close partialderiv;
use partialderiv;
read all var {df_da df_db df_dc df_dg df_dh df_di} into partialderiv;
close partialderiv;
Sigma = asycov;
rho_var = partialderiv*Sigma*t(partialderiv);
rho_SE = sqrt(rho_var);
l95=rhohat-1.96*rho_SE;
u95=rhohat+1.96*rho_SE;
p = (1-cdf("normal",abs(rhohat),0,rho_SE))*2;
NormalApproxOutput = j(1,5,.);
NormalApproxOutput[,1]=rhohat;
NormalApproxOutput[,2]=rho_SE;
NormalApproxOutput[,3]=l95;
NormalApproxOutput[,4]=u95;
NormalApproxOutput[,5]=p;
Outtitle={'Estimate' 'Std Error' '95% CI Lower Bound' '95% CI Upper
Bound' 'Pvalue'};
Varnames=t("&var1. and &var2.");
PRINT NormalApproxOutput[colname=Outtitle rowname=Varnames];
quit;

Ksharp · Posted 11-15-2020 06:20 AM

The IML is similar with data step.
The following code could translate into date step like:
proc iml;
use Covparms;
read all var {Estimate} into Cov_estimate;
close Covparms;
use Asycov;
read all var {CovP1 CovP2 CovP3 CovP4 CovP5 CovP6} into asycov;
close Asycov;
a = Cov_estimate[1];
b = Cov_estimate[2];
c = Cov_estimate[3];
g = Cov_estimate[4];
h = Cov_estimate[5];
i = Cov_estimate[6];
rhohat = (b+h)/sqrt((a+g)*(c+i));
df_da = -0.5*((b+h)*(c+i))/sqrt(((a+g)*(c+i))**3);
df_db = 1/sqrt((a+g)*(c+i));
df_dc = -0.5*(b+h)*(a+g)/sqrt(((a+g)*(c+i))**3);
df_dg = -0.5*(b+h)*(c+i)/sqrt(((a+g)*(c+i))**3);
df_dh = 1/sqrt((a+g)*(c+i));
df_di = -0.5*(b+h)*(a+g)/sqrt(((a+g)*(c+i))**3);

---------->>

data partialderiv ;
set Covparms end=last;
retain a b d g h i;
if _n_=1 then a=Estimate;
if _n_=2 then b=Estimate;
if _n_=3 then c=Estimate;
if _n_=4 then g=Estimate;
if _n_=5 then h=Estimate;
if _n_=6 then i=Estimate;

if last then do;
rhohat = (b+h)/sqrt((a+g)*(c+i));
df_da = -0.5*((b+h)*(c+i))/sqrt(((a+g)*(c+i))**3);
df_db = 1/sqrt((a+g)*(c+i));
df_dc = -0.5*(b+h)*(a+g)/sqrt(((a+g)*(c+i))**3);
df_dg = -0.5*(b+h)*(c+i)/sqrt(((a+g)*(c+i))**3);
df_dh = 1/sqrt((a+g)*(c+i));
df_di = -0.5*(b+h)*(a+g)/sqrt(((a+g)*(c+i))**3);

output;
end;
keep df_da df_db df_dc df_dg df_dh df_di;
run;

tka726 · Posted 11-30-2020 10:35 AM

Thanks so much!

Any idea how to code this part? I can not figure it out.

proc iml;
use Covparms;
read all var {Estimate} into Cov_estimate;
close Covparms;
use Asycov;
read all var {CovP1 CovP2 CovP3 CovP4 CovP5 CovP6} into asycov;
close Asycov;
a = Cov_estimate[1];
b = Cov_estimate[2];
c = Cov_estimate[3];
g = Cov_estimate[4];
h = Cov_estimate[5];
i = Cov_estimate[6];
rhohat = (b+h)/sqrt((a+g)*(c+i));
df_da = -0.5*((b+h)*(c+i))/sqrt(((a+g)*(c+i))**3);
df_db = 1/sqrt((a+g)*(c+i));
df_dc = -0.5*(b+h)*(a+g)/sqrt(((a+g)*(c+i))**3);
df_dg = -0.5*(b+h)*(c+i)/sqrt(((a+g)*(c+i))**3);
df_dh = 1/sqrt((a+g)*(c+i));
df_di = -0.5*(b+h)*(a+g)/sqrt(((a+g)*(c+i))**3);

create partialderiv var {df_da df_db df_dc df_dg df_dh df_di};
append;
close partialderiv;
use partialderiv;
read all var {df_da df_db df_dc df_dg df_dh df_di} into partialderiv;
close partialderiv;

Sigma = asycov;
rho_var = partialderiv*Sigma*t(partialderiv);
rho_SE = sqrt(rho_var);
l95=rhohat-1.96*rho_SE;
u95=rhohat+1.96*rho_SE;
p = (1-cdf("normal",abs(rhohat),0,rho_SE))*2;
NormalApproxOutput = j(1,5,.);
NormalApproxOutput[,1]=rhohat;
NormalApproxOutput[,2]=rho_SE;
NormalApproxOutput[,3]=l95;
NormalApproxOutput[,4]=u95;
NormalApproxOutput[,5]=p;
Outtitle={'Estimate' 'Std Error' '95% CI Lower Bound' '95% CI Upper
Bound' 'Pvalue'};
Varnames=t("&pltct. and &tegma.");
PRINT NormalApproxOutput[colname=Outtitle rowname=Varnames];
quit;

Ksharp · Posted 12-01-2020 07:44 AM

The obstacle is here:
rho_var = partialderiv*Sigma*t(partialderiv);

It is quadratic form .

data partialderiv ;
set Covparms end=last;
retain a b d g h i;
if _n_=1 then a=Estimate;
if _n_=2 then b=Estimate;
if _n_=3 then c=Estimate;
if _n_=4 then g=Estimate;
if _n_=5 then h=Estimate;
if _n_=6 then i=Estimate;

if last then do;
rhohat = (b+h)/sqrt((a+g)*(c+i));
df_da = -0.5*((b+h)*(c+i))/sqrt(((a+g)*(c+i))**3);
df_db = 1/sqrt((a+g)*(c+i));
df_dc = -0.5*(b+h)*(a+g)/sqrt(((a+g)*(c+i))**3);
df_dg = -0.5*(b+h)*(c+i)/sqrt(((a+g)*(c+i))**3);
df_dh = 1/sqrt((a+g)*(c+i));
df_di = -0.5*(b+h)*(a+g)/sqrt(((a+g)*(c+i))**3);

output;
end;
keep rhohat df_da df_db df_dc df_dg df_dh df_di; /*<--*/
run;
/***********************/
/***********************/
data _null_;
set partialderiv;
temp=catx(' ',df_da, df_db, df_dc, df_dg, df_dh, df_di);
call symputx('value',temp);
stop;
run;

data cov;
 set Asycov;
 array x{6} _temporary_ (&value);
 array c{*} CovP1-CovP6;
 i+1;
 do j=1 to dim(c);
  c{j}=c{j}*x{i}*x{j};
 end;
 drop i j;
run;

proc sql;
select sum(sum(CovP1, CovP2, CovP3 ,CovP4 ,CovP5, CovP6)) into : rho_var
 from cov;
quit;

data want;
set partialderiv;
rho_SE = sqrt(&rho_var);
l95=rhohat-1.96*rho_SE;
u95=rhohat+1.96*rho_SE;
p = (1-cdf("normal",abs(rhohat),0,rho_SE))*2;

label
rhohat='Estimate';
rho_SE='Std Error';
l95='95% CI Lower Bound';
u95='95% CI UpperBound';
p='Pvalue' ;
run;

Compute standard errors without using PROC IML

Re: Need help converting IML

Re: PROC IML Help!

Need help converting IML

Re: Need help converting IML

Re: Need help converting IML

Re: Need help converting IML