data oneway;
input trt y;
cards;
1  19.4
1  32.6
1  27.0
1  32.1
1  33.0
2  17.7
2  24.8
2  27.9
2  25.2
3  17.0
3  19.4
3   9.1
3  11.9
4  20.7
4  21.0
4  20.5
4  18.8
4  18.6
5  14.3
5  14.4
5  11.8
5  11.6
5  14.2
6  17.3
6  19.4
6  19.1
6  16.9
6  20.8
;
proc glm;
classes trt;
model y = trt/xpx i solution;
/* All possible pair-wise comparisons,
   even though excessive!
   Only 5 linearly independent, pre-planned
   comparisons should be made.
 N.B. in GLM the contrast statement pertains to the 
 specified factor, in this case trt, and the
 coefficients relate to trt  */
contrast ' trt 1 vs trt 2' trt 1 -1 0 0 0 0/e;
/* Note the option /e. This causes SAS GLM to print
out the whole k` matrix that it is using to construct
the requested contrast. Quite useful !  */
contrast ' trt 1 vs trt 3' trt 1 0 -1 0 0 0;
contrast ' trt 1 vs trt 4' trt 1 0 0 -1 0 0;
contrast ' trt 1 vs trt 5' trt 1 0 0 0 -1 0;
contrast ' trt 1 vs trt 6' trt 1 0 0 0 0 -1;
contrast ' trt 2 vs trt 3' trt 0 1 -1 0 0 0;
contrast ' trt 2 vs trt 4' trt 0 1 0 -1 0 0;
contrast ' trt 2 vs trt 5' trt 0 1 0 0 -1 0;
contrast ' trt 2 vs trt 6' trt 0 1 0 0 0 -1;
contrast ' trt 3 vs trt 4' trt 0 0 1 -1 0 0;
contrast ' trt 3 vs trt 5' trt 0 0 1 0 -1 0;
contrast ' trt 3 vs trt 6' trt 0 0 1 0 0 -1;
contrast ' trt 4 vs trt 5' trt 0 0 0 1 -1 0;
contrast ' trt 4 vs trt 6' trt 0 0 0 1 0 -1;
contrast ' trt 5 vs trt 6' trt 0 0 0 0 1 -1;
contrast ' SSR_m  ' trt 1 -1 0 0 0 0,
                    trt 1 0 -1 0 0 0,
                    trt 1 0 0 -1 0 0,
                    trt 1 0 0 0 -1 0,
                    trt 1 0 0 0 0 -1;

/* Estimates of differences */
estimate ' trt 1 - trt 2' trt 1 -1 0 0 0 0/e;
estimate ' trt 1 - trt 3' trt 1 0 -1 0 0 0;
estimate ' trt 1 - trt 4' trt 1 0 0 -1 0 0;
estimate ' trt 1 - trt 5' trt 1 0 0 0 -1 0;
estimate ' trt 1 - trt 6' trt 1 0 0 0 0 -1;
estimate ' trt 2 vs trt 3' trt 0 1 -1 0 0 0;
estimate ' trt 2 vs trt 4' trt 0 1 0 -1 0 0;
estimate ' trt 2 vs trt 5' trt 0 1 0 0 -1 0;
estimate ' trt 2 vs trt 6' trt 0 1 0 0 0 -1;
estimate ' trt 3 vs trt 4' trt 0 0 1 -1 0 0;
estimate ' trt 3 vs trt 5' trt 0 0 1 0 -1 0;
estimate ' trt 3 vs trt 6' trt 0 0 1 0 0 -1;
estimate ' trt 4 vs trt 5' trt 0 0 0 1 -1 0;
estimate ' trt 4 vs trt 6' trt 0 0 0 1 0 -1;
estimate ' trt 5 vs trt 6' trt 0 0 0 0 1 -1;
/* Estimates of fitted values */
estimate ' mean + trt 1 ' intercept 1 trt 1 0 0 0 0 0;
estimate ' mean + trt 2 ' intercept 1 trt 0 1 0 0 0 0;
estimate ' mean + trt 3 ' intercept 1 trt 0 0 1 0 0 0;
estimate ' mean + trt 4 ' intercept 1 trt 0 0 0 1 0 0;
estimate ' mean + trt 5 ' intercept 1 trt 0 0 0 0 1 0;
estimate ' mean + trt 6 ' intercept 1 trt 0 0 0 0 0 1;
/* Least Squares Means, standard errors and probability
   of differences between pairs of treatments */
lsmeans trt/stderr pdiff;
run;

proc iml; /* Explicit construction of matrices, etc, using PROC IML,
         IML = Interactive Matrix Language */
reset print;

x = {1 1 0 0 0 0 0,
     1 1 0 0 0 0 0,
     1 1 0 0 0 0 0,
     1 1 0 0 0 0 0,
     1 1 0 0 0 0 0,
     1 0 1 0 0 0 0,
     1 0 1 0 0 0 0,
     1 0 1 0 0 0 0,
     1 0 1 0 0 0 0,
     1 0 1 0 0 0 0,
     1 0 0 1 0 0 0,
     1 0 0 1 0 0 0,
     1 0 0 1 0 0 0,
     1 0 0 1 0 0 0,
     1 0 0 1 0 0 0,
     1 0 0 0 1 0 0,
     1 0 0 0 1 0 0,
     1 0 0 0 1 0 0,
     1 0 0 0 1 0 0,
     1 0 0 0 1 0 0,
     1 0 0 0 0 1 0,
     1 0 0 0 0 1 0,
     1 0 0 0 0 1 0,
     1 0 0 0 0 1 0,
     1 0 0 0 0 1 0,
     1 0 0 0 0 0 1,
     1 0 0 0 0 0 1,
     1 0 0 0 0 0 1,
     1 0 0 0 0 0 1,
     1 0 0 0 0 0 1};

y = {19.4,
     32.6,
     27.0,
     32.1,
     33.0,
     17.7,
     24.8,
     27.9,
     25.2,
     24.3,
     17.0,
     19.4,
      9.1,
     11.9,
     15.8,
     20.7,
     21.0,
     20.5,
     18.8,
     18.6,
     14.3,
     14.4,
     11.8,
     11.6,
     14.2,
     17.3,
     19.4,
     19.1,
     16.9,
     20.8};

use oneway;     /* Another way to get the data in IML, from
                   a SAS data set */
read all var{trt y} into plants;

x1 = plants[,1];
nobs = nrow(x1);
x0 = j(nobs,1,1);
x = x0 || design(x1);
y = plants[,2];


xtx = x` * x; /* X'X */
xty = x` * y; /* X'Y */
invxtx = ginv(xtx); /* generalised inverse of X'X */
b = invxtx * xty; /* solution vector */
tss = y` * y; /* Total Sums of Squares */
sumy = sum(y); /* sum y values */
nobs = nrow(x); /* number of rows in X, = N = number of observations */
ssr = b` * xty; /* SSR, Sums of Squares for the Model */
ybar = sumy/nobs; /* average of the y's */
cf = nobs * ybar * ybar; /* Corection Factor for the mean */
ssrm = ssr - cf; /* SSR_m */
sse = tss - ssr; /* Residual Sums of Squares */
dfe = nobs - ncol(x) + 1; /* Residual Degrees of Freedom */
mse = sse/dfe; /* Mean Square Error */

/* Estimates of fitted values and their standard errors */

k1 = {1, 1, 0, 0, 0, 0, 0};
kb = k1` * b;
kgk = k1` * invxtx * k1;
sv = kgk * mse;
se = sqrt(sv);
k2 = {1, 0, 1, 0, 0, 0, 0};
kb = k2` * b;
kgk = k2` * invxtx * k2;
sv = kgk * mse;
se = sqrt(sv);
k3 = {1, 0, 0, 1, 0, 0, 0};
kb = k3` * b;
kgk = k3` * invxtx * k3;
sv = kgk * mse;
se = sqrt(sv);
k4 = {1, 0, 0, 0, 1, 0, 0};
kb = k4` * b;
kgk = k4` * invxtx * k4;
sv = kgk * mse;
se = sqrt(sv);
k5 = {1, 0, 0, 0, 0, 1, 0};
kb = k5` * b;
kgk = k5` * invxtx * k5;
sv = kgk * mse;
se = sqrt(sv);
k6 = {1, 0, 0, 0, 0, 0, 1};
kb = k6` * b;
kgk = k6` * invxtx * k6;
sv = kgk * mse;
se = sqrt(sv);

/* Estimates of differences between treatments and the standard */
/* errors of the differences */

k12 = {0, 1, -1, 0, 0, 0, 0};
kb = k12` * b;
kgk = k12` * invxtx * k12;
sv = kgk * mse;
se = sqrt(sv);
t12 = kb/se;
pr = 2*(1 - probt(abs(t12),dfe));
invkk = inv(kgk);
ss12 = kb` * invkk * kb;
df = ncol(k12);
ms12 = ss12/df;
f = ms12/mse;
pr = 1 - probf(f,df,dfe);
ss = kb` * inv(kgk) * kb;
ms = ss/1;
f12 = ms/mse;

k13 = {0, 1, 0, -1, 0, 0, 0};
kb = k13` * b;
kgk = k13` * invxtx * k13;
sv = kgk * mse;
se = sqrt(sv);
t13 = kb/se;
ss = kb` * inv(kgk) * kb;
ms = ss/1;
f13 = ms/mse;

k14 = {0, 1, 0, 0, -1, 0, 0};
kb = k14` * b;
kgk = k14` * invxtx * k14;
sv = kgk * mse;
se = sqrt(sv);
t14 = kb/se;
ss = kb` * inv(kgk) * kb;
ms = ss/1;
f14 = ms/mse;

k15 = {0, 1, 0, 0, 0, -1, 0};
kb = k15` * b;
kgk = k15` * invxtx * k15;
sv = kgk * mse;
se = sqrt(sv);
t15 = kb/se;
ss = kb` * inv(kgk) * kb;
ms = ss/1;
f15 = ms/mse;

k16 = {0, 1, 0, 0, 0, 0, -1};
kb = k16` * b;
kgk = k16` * invxtx * k16;
sv = kgk * mse;
se = sqrt(sv);
t16 = kb/se;
ss = kb` * inv(kgk) * kb;
ms = ss/1;
f16 = ms/mse;

kp = {0  1  -1  0  0  0  0,
      0  1  0  -1  0  0  0,
      0  1  0  0  -1  0  0,
      0  1  0  0  0  -1  0,
      0  1  0  0  0  0  -1};
kb = kp * b;
kgk = kp * invxtx * kp`;
ss = kb` * inv(kgk) * kb;
dftrt = nrow(kp);
ms = ss/dftrt;
f = ms/mse;
pr = 1 - probf(f,dftrt,dfe);
quit;