The MxMult macro always attempts to check for a later version of itself. If it is unable to do this (such as if there is no active internet connection available), the macro issues the following message in the log:

   NOTE: Unable to check for newer version of the MXMULT macro.

The computations performed by the macro are not affected by the appearance of this message. However, this check can be avoided by specifying the following statement in your SAS session. This can be useful if your machine has no connection to the internet.

   %let newchk=x;

Download and save mxmult.sas

   %inc "<location of your file containing the MxMult macro>";

You can now call the MxMult macro during your SAS session.

The MxMult macro requires that you specify one or two data set names in the first two arguments. Each data set should contain numeric variables that hold the matrix to be operated on. In the following, the data sets in the first two arguments are referred to as matrix A and matrix B.

To transpose, diagonalize, or invert a single matrix, specify only one data set name. To do any of these operations on each of two matrices and then multiply, add, or subtract them, specify two data set names. Some examples:

%mxmult(a,diag=y)

Diagonalizes the 1 row (or column) vector in data set a (matrix A) to create a square matrix with the vector on the main diagonal.

%mxmult(a,b,t=y n)

Transposes the matrix in data set a (matrix A) and multiplies it by the matrix in data set b (matrix B): A'B

%mxmult(x,y,inv=y n,t=n y)

Inverts the square, nonsingular matrix in data set x (matrix A) and multiplies it by the transpose of the matrix in data set y (matrix B): X^-1Y'

%mxmult(x,y,op=a)

Adds the matrix in data set x (matrix A) and the matrix in data set y (matrix B): A+B

Optional parameters

Only one operation (transpose, diagonal, invert) can be performed on each input matrix.

t=N|Y <N|Y>

Requests that one or both matrices be transposed before multiplication. By default, t=N N indicating that neither is transposed. The first N or Y refers to matrix A. The second N or Y (if two matrices are specified) refers to matrix B. If you specify t= for a matrix, you cannot also specify diag=, vdiag=, or inv= for that matrix.

diag=N|Y <N|Y>

Requests that one or both matrices be diagonalized before multiplication. By default, diag=N N indicating that neither is diagonalized. The first N or Y refers to matrix A. The second N or Y (if two matrices are specified) refers to matrix B. If you specify diag= for a matrix, you cannot also specify t=, vdiag= or inv= for that matrix.

vdiag=N|Y <N|Y>

For one or both matrices, Y requests extraction of the main diagonal to a column vector. The matrix must be square. By default, vdiag=N N indicating no diagonal extraction for either matrix. The first N or Y refers to matrix A. The second N or Y (if two matrices are specified) refers to matrix B. If you specify vdiag= for a matrix, you cannot also specify t=, diag= or inv= for that matrix.

inv=N|Y <N|Y>

Requests that one or both matrices be inverted before multiplication. By default, inv=N N indicating that neither is inverted. The first N or Y refers to matrix A. The second N or Y (if two matrices are specified) refers to matrix B. If you specify inv= for a matrix, you cannot also specify t=, diag=, or vdiag= for that matrix.

op=A|S|M

When both A and B matrices are specified, op= indicates whether they should be added, subtracted, or multiplied. When only one matrix is specified, op= is ignored. The default is multiplication, op=M. Matrices are checked for compatibility. See Details for more on compatibility.

out=data-set-name

Specifies the name of the data set created containing the results of the requested operations. By default, out=AxB.

prefix=string

Specifies a prefix for the variable names of the matrix columns in the out= data set. When the result is scalar, it also provides the name of a macro variable (with 1 appended) containing the scalar result. See Details for more information. By default, prefix=Prod.

You can use the MxMult macro to perform certain operations on one or two matrices. The columns of each matrix should be contained in numeric variables in a data set. Specify one or two data set names in the first two arguments of the macro. For example, if two matrices are contained in the numeric variables of data sets A and B, and assuming that the matrices are compatible, then %mxmult(A,B) creates the product, A × B, and stores the resulting matrix in data AxB with column names Prod1, Prod2, ... .

Prior to multiplying, adding, or subtracting matrices A and B, each matrix can optionally be transposed, diagonalized, inverted, or have its main diagonal saved as a vector. Only one of these operations can be done on each matrix. In order to diagonalize a matrix, the matrix must be a one row or one column matrix (that is, a vector). If the vector has p elements, the resulting diagonalized matrix has p rows and p columns with the p elements of the vector on the main diagonal. In order to invert a matrix or extract its main diagonal, it must be a square matrix (same number of rows and columns). For inversion, it must not be singular. If diagonal extraction is requested, the main diagonal is saved in a column vector. In order for two matrices to be multiplied, added, or subtracted, they must be compatible following any transposition, diagonalization, inversion, or diagonal extraction. To be compatible for multiplication, the number of columns of matrix A must equal the number of rows of matrix B. In order to be compatible for addition or subtraction, the matrices must have the same numbers of rows and columns. A Note appears in the SAS Log giving the final dimensions of each matrix.

Neither matrix can be a scalar (one row and one column). Prior to performing any operation on either matrix, the data sets are cleaned to remove any character variables and any observations that contain a missing value. Note that observations with missing values in character variables will cause the deletion of those observations.

Saved results from the MxMult macro

The MxMult macro produces no displayed results. The final computed matrix is saved in the out= data set, named AxB by default. The variable names in the out= data set all begin with the string specified in prefix= and are appended with 1, 2, 3, and so on. By default, prefix=Prod, so the variable names are Prod1, Prod2, Prod3, and so on.

If the final matrix is a scalar, then it is also saved in a macro variable available after the macro terminates. The name of the macro variable is the same as the single variable in the out= data set as described above. That is, the name begins with the string specified in prefix= and is appended with 1. By default, prefix=Prod, so a scalar result is saved in macro variable Prod1. When a macro variable is created, a Note appears in the SAS Log giving its name.

Example 1: Diagonalize and invert a matrix

You can invert any square, nonsingular matrix. The data step below first creates a 1×4 row vector and stores the values in four numeric variables in data set MX14.

data MX14; 
   array x (4) (1:4);
   run;
proc print noobs;
   run;

The first macro call creates a square diagonal matrix from the row vector and saves it in data set Diag. Note that a column vector could also be used.

%mxmult(MX14,diag=y,out=Diag)
proc print data=Diag noobs;
   run;

The second macro call inverts the Diag matrix. Since it is a diagonal matrix, its inverse is simply a diagonal matrix with the reciprocals of the values on the diagonal.

%mxmult(Diag,inv=y,out=Inv)
proc print data=Inv noobs;
   format _all_ 4.2;
   run;

This final macro call uses vdiag= to extract the main diagonal of the square matrix above and saves it as a column vector in out= data set Vdiag.

%mxmult(Inv,vdiag=y,out=Vdiag)
proc print data=Vdiag noobs;
   format _all_ 4.2;
   run;

Example 2: Estimate regression parameters

You can compute more complex matrix expressions by combining multiple calls of the macro. To illustrate, the following computes the estimates of a linear regression model. This example estimates the parameters of the model

E(weight) = b0 + b1*height + b2*age

using the data in SASHELP.CLASS.

The matrix formula for the parameter estimates is

inv(X'X)X'Y .

To begin, the columns of the X matrix are stored in data set X. Note that a column of 1 values is added for the intercept. The LENGTH statement orders the columns. Similarly, the response values are stored in data set Y.

data X(keep=one height age); 
   length one height age 8.;
   set sashelp.class;
   one=1;
   run;
data Y(keep=weight);
   set sashelp.class;
   run;

To compute the estimating formula, it is broken into a series of multiplications of two matrices at a time using any needed operation on each matrix. The first macro call computes X'X and stores it in data set XPX. The next call inverts X'X, multiplies it by X', and stores the result in data set XPXinvX. The formula is completed by multiplying the last result by Y. The final estimates are stored in data set BetaHat.

%mxmult(X,X,t=y n,out=XPX)
%mxmult(XPX,x,inv=y n,t=n y,out=XPXinvX)
%mxmult(XPXinvX,y,out=BetaHat)

These statements display the estimates computed using the macro and compares them to the estimates computed using PROC GLM.

proc print data=BetaHat noobs; 
   title "Parameter estimates";
   run;
proc glm data=sashelp.class;
   model weight=height age;
   output out=glmout p=pred;
   ods select parameterestimates;
   quit;

Example 3: Compute predicted values

Using the data and parameter estimates in the previous example, the following macro call computes predicted values as Xβ. The predicted values are saved in variable Prod1 in data set AxB. The DATA step combines them with the predicted values computed above by PROC GLM and displays them both showing their agreement. The first five observations are shown.

%mxmult(x,betahat)
data result;
   merge glmout axb;
   keep pred prod1;
   run;
proc print data=result(obs=5) noobs;
   title "Predicted values from GLM and MxMult";
   run;

Example 4: Compute Wald chi-square for a contrast of model parameters

Given a column vector of model parameters, β, the covariance matrix of model parameters, Σ, and one or more row vectors in a matrix, L, a Wald chi-square statistic testing the hypotheses Lβ=0 is

(Lβ)' (LΣL')^-1 Lβ .

The following statements fit the model appearing in the Poisson regression example in the Getting Started section of the GENMOD procedure documentation. The parameter estimates and covariance matrix are saved in data sets using an ODS OUTPUT statement. The CONTRAST statement tests that the parameters of the model predictors are simultaneously zero, providing an overall test of the model. The DATA L step creates a hypothesis matrix in which each row, when multiplied by β, selects one predictor parameter. As a result, Lβ represents the hypothesis that all predictor parameters are zero.

proc genmod data=insure;
   class car age / param=ref;
   model c = car age / dist=poisson offset=ln covb;
   contrast 'model' car 1 0, car 0 1, age 1 / wald;
   ods output covb=Cov
     parameterestimates=Beta(where=(parameter ne 'Scale')); 
   run;
data Beta(keep=estimate); 
   set Beta; 
   run;
data L;
   input L1-L4;
   datalines;
0 1 0 0
0 0 1 0
0 0 0 1
;

This table from GENMOD shows the test performed by CONTRAST statement.

These macro calls reproduce the Wald chi-square statistic computed by the CONTRAST statement above.

/* L*Beta */
%mxmult(L,Beta,out=LBeta)
/* L*Cov */
%mxmult(L,Cov,out=LCov)
/* L*Cov*L' */
%mxmult(LCov,L,t=n y,out=LCovL)
/* (L*Beta)'*inv(L*Cov*L') */
%mxmult(LBeta,LCovL,t=y n,inv=n y,out=LBLCLinv)
/* (L*Beta)'*inv(L*Cov*L')*(L*Beta) */
%mxmult(LBLCLinv,LBeta,out=CSQ,prefix=Chisq)

These statements display the resulting chi-square statistic and its p-value. The statistic has 3 degrees of freedom due to the three independent rows in L.

data chisq; set CSQ;
   prob=sdf('CHISQ',Chisq1,3);
   run;
proc print noobs;
   format prob pvalue6.;
   title "Wald chi-square for model test";
   run;

The results of the matrix operations below agree with the test from the CONTRAST statement above.

Example 5: Confidence interval for ratio of random variables using Fieller's theorem

SAS Note 56476 presents a SAS/IML-based macro to compute a confidence interval for a ratio of random variables using Fieller's theorem. The same can be done without SAS/IML by using the MxMult macro. See the SAS Note for details on the computation and an example. The following reproduces the example results for the first ratio using MxMult macro calls and a DATA step.

Since each of the matrix multiplications shown below produce a scalar result, you can take advantage of the macro saving each result in a macro variable. The name of the macro variable is the prefix= string appended with 1. Notice the references to the macro variables in the final DATA step that computes the confidence limits.

proc logistic data=assay;
   class drug / param=glm;
   model ndead/total = drug ldose / noint covb;
   ods output covb=Cov
       parameterestimates=Beta(keep=estimate);
   run;
data k;
   array k (3) (-1 0 0);
   run;
data h;
   array h (3) (0 0 1);
   run;
%mxmult(k,beta,prefix=kb)
%mxmult(h,beta,prefix=hb)
%mxmult(h,cov,out=hcov)
%mxmult(k,cov,out=kcov)
%mxmult(hcov,h,t=n y,prefix=hch)
%mxmult(kcov,h,t=n y,prefix=kch)
%mxmult(kcov,k,t=n y,prefix=kck)
data abc;
   tsq=probit(1-0.05/2)**2;
   ratio=&kb1/&hb1;
   a=&hb1**2-tsq*&hch1;
   b=2*(tsq*&kch1-&kb1*&hb1);
   c=&kb1**2-tsq*&kck1;
   disc=((b**2)-4*a*c);
   if (disc<=0 | a<=0) then
     put "ERROR: Confidence interval can't be computed";
   else do;
     sroot=sqrt(disc);
     l_b=((-b)-sroot)/(2*a);
     u_b=((-b)+sroot)/(2*a);
   end;
   run;
proc print noobs;
   var ratio l_b u_b;
   run;

The results from the above statements replicate the confidence interval for the first ratio shown in the SAS Note.