Rounding issue

[mgx]

Dear SAS community,

I am working with an imputed dataset and first performing a sandwich estimator to compute the within-imputation variances on each imputation before combining the results. My proc iml code for the sandwich estimator looks as follows:

%macro sandwich_imput(i=);
data hessian;
set mvp.hessian_imput;
where _imputation_=&i;
drop row Parameter _imputation_;
run;
data gradient;
set mvp.grad_imput;
where _imputation_=&i;
drop observation _imputation_;
run;
proc iml;
use hessian; read all into J; close hessian;
use gradient; read all into G; close gradient;
use parms; read all into P; close parms;
K=t(G)*(G);
Sigma=round(inv(J)*K*inv(J),0.001);
create Sigma from Sigma;append from Sigma;
quit;
data cov;
set Sigma cov;
run;
%mend sandwich_imput;

After adding the appropriate extra variables, I use the cov dataset in PROC MIANALYZE. However, due to rounding errors I get the following error: 

ERROR: Within-imputation COVB matrix is not symmetric for _Imputation_= 1 in the input COVB= data set.

The rounded variance-covariance matrix is symmetric, but it seems that SAS internally stores the non-rounded matrix and uses that matrix in PROC MIANALYZE. Is there any way to circumvent this problem?

[FreelanceReinh]

Hello @mgx,

---

@mgx wrote:

ERROR: Within-imputation COVB matrix is not symmetric for _Imputation_= 1 in the input COVB= data set.

The rounded variance-covariance matrix is symmetric, but it seems that SAS internally stores the non-rounded matrix and uses that matrix in PROC MIANALYZE. Is there any way to circumvent this problem?

---

SAS will use the matrix given in the "input COVB= data set" (which you named COV_IMPUT in your other thread). So I would look at this dataset, e.g., with PROC PRINT:

proc print data=cov_imput;
where _imputation_=1;
*format _numeric_ best16.;
run;

to see where (and by how much) the matrix deviates from symmetry. Use the FORMAT statement (commented out above) if you see matrix elements with more than three decimals. Actually this should not happen because, as far as I see, all these numbers arise from

Sigma=round(inv(J)*K*inv(J),0.001);

So if a matrix in COV_IMPUT deviates from symmetry, then the original Sigma must be asymmetric as well, unless you introduced the asymmetry when creating COV_IMPUT from COV. Examining that original Sigma or the code creating COV_IMPUT would then be the next step in the investigation.

[mgx]

Hi FreelanceReinh, 

My matrix is symmetric until about 9 decimal points or so. Hence, it is a rounding issue and I am looking for a way to decrease the number of decimal points of the matrix.

[FreelanceReinh]

---

@mgx wrote:

Hi FreelanceReinh, 

My matrix is symmetric until about 9 decimal points or so. Hence, it is a rounding issue and I am looking for a way to decrease the number of decimal points of the matrix.

---

I assumed that

Sigma=round(inv(J)*K*inv(J),0.001);

ensures that all values are rounded to three decimal places (couldn't test it, though, for lack of a SAS/IML license).

Can you show one of the pairs of off-diagonal elements in COV_IMPUT that make the variance-covariance matrix for _imputation_=1 asymmetric?

[Rick_SAS · Accepted solution]

I don't know if the issue is from the Hessian matrix that you are reading in or from the rounding that you are manually performing in your program.

Why do you think you need to round the Sigma matrix to 0.001?  If any element is less than 0.0005 in absolute value, it will be truncated to 0. Is that what you want?

Here are some IML helper functions that you can use to diagnose and fix the issue:

/* test symmetry */
start TestSymmetry(M);
   diff = max(abs(M - M`));
   name = strip(parentname("M"));
   labl = "Max Abs Diff: " + name + " - T(" + name + ")";
   print diff[L=labl];
finish;

/* make a matrix symmetric */
start MakeSymmetric(M);
   return (M + M`) / 2;
finish;

The first module prints the maximum absolute difference of any unsymmetric terms. The second module returns a symmetric matrix that is close to the input matrix. 

I suggest you test the input Hessian (J) for symmetry by using RUN TestSymmetry(J);  If the module prints 0, then the input is not a problem/

My suspicion is that the problem in in the computation of Sigma. If J is an ill-conditioned matrix, then computing the inverse is an unstable operation and Sigma might not be symmetric. I suggest you try

S = inv(J)*K*inv(J);
run TestSymmetry(S);
Sigma= MakeSymmetric(S);  /* ROUND here, if required */

I didn't use the ROUND function because I don't understand why it is needed. If required, you can add 
Sigma = round(Sigma, 0.001);
after the call to the MakeSymmetric function.