Let's say that I have a numeric data matrix with columns `w, x, y, z` and I also want to add in the columns that are equivalent to `w*x, w*y, w*z, x*y, x*z, y*z` since I want my covariate matrix to include all pairwise interactions.

Is there a clean and effective way to do this?