*Fit a smoothing spline using PROC MIXED; 

options ls=72;
data diabetes;
  infile 'diabetes.dat';
  input age basedef cpep;
lncpep=log(cpep);


* Identify distinct values of age;
data knots;
 set diabetes;

proc sort data=knots;
by age;

data knots;
 set knots (keep=age);
 by age;
 if first.age;


proc iml;

use diabetes;
read all into w;

use knots;
read all into t0;

t       = w[,1];
y       = w[,4];
n       = nrow(w);
nk      = nrow(t0);
h       = j(nk-1,1,0);
q       = j(nk,nk-2,0);
rr      = j(nk-2,nk-2,0);
g       = j(nk-2,nk-2,0);
l       = j(nk,nk-2,0);
b       = j(nk,nk-2,0);
z       = j(n,nk-2,0);

*calculate h;
do i=1 to nk-1;
 h[i]=t0[i+1]-t0[i];
end;

*calculate Q;
do j=1 to nk-2;
 q[j,j]=1/h[j];
 q[j+1,j]=-1/h[j+1]-1/h[j];
 q[j+2,j]=1/h[j+1];
end;

*calculate R;
do j=1 to nk-2;
 rr[j,j]=(h[j]+h[j+1])/3;
end;
rr[1,2]=h[2]/6;
do j=2 to nk-3;
 rr[j,j-1]=h[j]/6;
 rr[j,j+1]=h[j+1]/6;
end;
rr[nk-2,nk-3]=h[nk-2]/6;

*calculate G;
g=t(half(rr));

*calculate L;
l=q*t(inv(g));

*calculate B;
b=l*inv(t(l)*l);

*find index of element e in vector t0, i.e, calculate the incidence matrix N;
start index(e,t0);
 do i=1 to nrow(t0);
  if e=t0[i] then j=i;
 end;
 return(j);
finish;

*define design matrix Z=NB for introduced random effects;
do i=1 to n;
 do j=1 to nk-2;
  z[i,j]=b[index(t[i],t0),j];
 end;
end;

*create SAS data set for PROC MIXED;
datmtr=y||t||z;
create diabetes2 from datmtr;
append from datmtr;

quit;

data diabetes2;
 set diabetes2;
 rename col1=lncpep col2=age col3-col37=z1-z35;

proc mixed data=diabetes2 covtest;
 model lncpep=age/s outpred=fit;
 random z1-z35/type=toep(1) s;


proc sort data=fit;
by age;

data fit;
 set fit (keep=age Pred StdErrPred);
 by age;
 if first.age;
file "spline_mixed.out";
put age Pred StdErrPred;

proc print data=fit;
var age Pred StdErrPred;

*Fit a smoothing spline via mixed models using macro spmm.mac;
%include '$HOME/gamm/Daowen_program/spmm.mac';

%spmm(data=diabetes,dep=lncpep,smthvar=age,outband=ci, print=y);


data one;
set ci;
file 'spmm_fit.out';
put age lncpep F_Se F_L95 F_U95  B_Se  B_L95 B_U95;