options ls=80 nodate pageno=1 NOFMTERR; 
 
 *  Class13_11.sas  segmented regression   20 Oct 2011;

 Proc SGplot data=ph6470.segmented_reg;
   scatter x=X y=Y;
 Proc SGplot data=ph6470.segmented_reg;
   loess x=X y=Y / LINEATTRS = ( color = blue);

 *** assume breakpoint = 60 (true = 55);

  Proc NLin data=ph6470.segmented_reg;
    parms b1=10.0 b2=0.5  b3=0  ;
    if (X < 60.0) then mean=b1 + b2*X;
    else mean= b1 + (b2 - b3)*60.0 + b3*X;
    model Y = mean;
    output out=A predicted=yhat;

  Proc SGplot data=A;
    scatter x=X y=Y;
    series x=X y=yhat / LINEATTRS = ( color = red thickness=3);


  proc sort data=A;
    by X;

  Proc SGplot data=A;
   scatter x=X y=Y;
    series x=X y=yhat;

  Proc SGplot data=A;
   scatter x=X y=Y;
    series x=X y=yhat / LINEATTRS = ( color = red thickness=3);
   loess x=X y=Y / LINEATTRS = ( color = blue);

 run; quit;

  *** estimate breakpoint d (true = 55);

  Proc NLin data=ph6470.segmented_reg;
    parms b1=10.0 b2=0.5  b3=0 d=60.0 ;
    if (X < d) then mean=b1 + b2*X;
    else mean= b1 + (b2 - b3)*d + b3*X;
    model Y = mean;
    output out=B predicted=yhat ;

  proc sort data=B;
    by X;

  Proc SGplot data=B;
   scatter x=X y=Y;
    series x=X y=yhat / LINEATTRS = ( color = red thickness=3);

  Proc SGplot data=B;
   scatter x=X y=Y;
    series x=X y=yhat / LINEATTRS = ( color = red thickness=3);
   loess x=X y=Y / LINEATTRS = ( color = blue);


 run; quit;