data metab;
input genotype $ nsdl metab;
cards;
A 9 0.0231
B 6 0.0275
B 7 0.0521
A 9 0.0001
A 6 0.0063
B 6 0.0138
B 8 0.1061
A 5 0.0482
;
/* ols estimates */
proc glm;
class genotype;
model metab = genotype;
lsmeans genotype /stderr;
estimate 'diff B-A' genotype -1 1;
run;
/* weighted ls is easy */
/* the variable in the weight statement is W_i */
/* For Aitken with diagonal V, W_i is 1/sigma_i*/
proc glm;
class genotype;
model metab = genotype;
weight nsdl;
lsmeans genotype /stderr;
estimate 'diff B-A' genotype -1 1;
run;
/* for more general V matrices, can sometimes force SAS procedures to give you */
/* what you want. In general, need to work with IML */
proc iml;
/* read the metab data in */
use metab;
read all var {metab} into y;
read all var {genotype} into genotype;
read all var {nsdl} into nsdl;
nobs = nrow(genotype);
/* create full rank X matrix, using set first to 0 */
d = design(genotype);
x = j(nobs,1,1) || d[,2];
/* now create V matrix */
v = diag(1/nsdl);
vi = inv(v);
/* and C matrix */
c = {1 0, 0 1};
/* ready to compute what we need */
bh = ginv(x` * vi * x) * x` * vi * y;
s2 = (y - X * bh)` * vi * (y - X*bh) / (nobs - 2);
s = sqrt(s2);
cb = c * bh;
vcb = s2 * c * ginv(x`*vi*x) * c`;
se = sqrt(vcb[1,1] || vcb[2,2]);
print(cb);
print(vcb);
print(se);
print(s);
/* or create z and w and save for glm */
call eigen(eval, evect, V);
vis = evect * diag(1/sqrt(eval)) * evect`;
print(v);
print(vis);
w = vis * x;
z = vis * y;
/* to create a SAS data set from IML, need to collect everything into 1 matrix */
wz = w || z;
create newmetab from wz; /* create the data set */
append from wz; /* and append all obs */
quit;
proc print data = newmetab; /* contents of newmetab */
run;
proc glm;
model col3 = col1 col2 / noint; /* col3 is the response (z), cols 1 and 2 are W */
/* /noint suppresses the default intercept */
run;
/* possibly useful ideas to construct other forms of V matrices */
proc iml;
/* for an ar structure matrix, use toeplitz of powers of rho */
rho = 0.7;
power = 0:7; /* vector of 0, 1, 2, ... 7 */
rhop = rho ## power; /* vector of 1, rho, rho^2, ... rho^6 */
v = toeplitz(rhop); /* banded matrix based on rhop */
print(v);
/* for a compound symmetric structure, I don't have any great ideas */
/* if equal numbers of obs per group, can use kronecker products */
/* example is for 5 groups of 2 obs each */
/* vc matrix is block diagonal with 5 blocks, each [1+k 1, 1, 1+k] */
k = {2,2}; /* ratio of s^2_p / s^2_c */
vone = j(2,2,1) +diag(k); /* one block of the vc matrix */
print(vone);
v = i(5) @ vone ; /* @ is the kronecker product operator */
/* block diagonal matrix is identity matrix @ vone */
print(v);
quit;