/* "Randomized Dose-Response Study of Rosuvastin in Japanese Patients with 
   Hypercholesterolemia" Saito, et al, Journal of Atherosclerosis, 2003, 10:329-336

   Weighted Least Squares:
   Y_i=%Change in LDL Cholesterol at 12 weeks for Dose_i
   X_i=LN(Dose_i)
   W_i=SQRT(r_i)
   r_i=Number of subjects at Dose_i
*/


options nodate nonumber ps=54 ls=80;
goptions reset=all dev=pdf gsfname=output gsfmode=replace;
data one;
input i 1 dLDL 4-8  DOSE 11-13 r 16-17;
lnDOSE=log(DOSE);
cards;
1  -35.8    1  15
2  -45.0  2.5  17
3  -52.7    5  12
4  -49.7   10  14
5  -58.2   20  18
6  -66.0   40  13
;
run;

proc reg;
weight r;
model dLDL=lnDOSE / p r;
output out=regout p=yhat;
run;


filename output 'wlsdr.pdf';

symbol1 c=black l=1 i=join;
symbol2 c=black v=:;

proc gplot;
plot yhat*lnDOSE=1 dLDL*lnDOSE=2 / overlay frame;
title 'Weighted Least Squares -- Dose-Response Cholesterol Study';
run;

quit;

proc iml;
y={-35.8,-45.0,-52.7,-49.7,-58.2,-66.0};
x0=j(6,1,1);
x1r={1, 2.5, 5, 10, 20, 40};
w2={15,17,12,14,18,13};
x1=log(x1r);
x=x0||x1;
w=sqrt(w2);
w=diag(w);

xs=w*x;
ys=w*y;

betahatw=inv(xs`*xs)*xs`*ys;

es=ys-xs*betahatw;
ssew=es`*es;
msew=ssew/(6-2);

vbetahatw=diag(msew*inv(xs`*xs));

sebetahatw=sqrt(vbetahatw)*j(2,1,1);


yhatw=x*betahatw;

ew=y-yhatw;

print betahatw sebetahatw yhatw ew;

run;

stop;