topic Re: How to do Regression on Broken/Turning Line?! in SAS Programming

How to do Regression on Broken/Turning Line?!

hellohere — Wed, 07 Sep 2022 20:06:36 GMT

%let slope1=-2; %let slope2=5; %let turner=30;

data indata;
do i=1 to 100;
	if i<&turner. then y=5*ranuni(i)+i*&slope1.;
	else y=10*ranuni(i)+i*&slope2.-&turner.*&slope2.+&turner.*&slope1.;
	output;
end;
run;quit;

ods listing gpath="%sysfunc(getoption(work))";
proc sgplot data=indata;
scatter x=i y=y; 
run;quit;

I have observations on a broken line alike the one above.

How to do a fitting to find the parameters (slope1, slope2 and turner)?! Thanks,

Re: How to do Regression on Broken/Turning Line?!

hellohere — Wed, 07 Sep 2022 20:15:07 GMT

I mean to fit stepwise line once for all, since turner[the turning point] is uncertain in general.

Re: How to do Regression on Broken/Turning Line?!

PaigeMiller — Wed, 07 Sep 2022 20:57:06 GMT

Here's an example with a curved portion and a linear portion; just replace the curved portion with a straight line

https://documentation.sas.com/doc/en/pgmsascdc/9.4_3.4/statug/statug_nlin_examples01.htm

Re: How to do Regression on Broken/Turning Line?!

hellohere — Thu, 08 Sep 2022 01:00:41 GMT

proc nlin data=indata;
   parms slope1=0 slope2=0 x0=50;

   if (x < x0) then
        mean = x*slope1;
   else mean = x*slope2-x0*slope2+x0*slope1;
   model y = mean;

   output out=indata_out predicted=yp;
run;
proc sgplot data=indata_out;
scatter x=x y=y; 
scatter x=x y=yp; 
run;quit;

Ye, it works. But is set the starting as slope1=-0.1 slope2=0.1, the outcomes are quite off.

Re: How to do Regression on Broken/Turning Line?!

PaigeMiller — Thu, 08 Sep 2022 10:12:57 GMT

Try other starting values

Re: How to do Regression on Broken/Turning Line?!

hellohere — Fri, 16 Sep 2022 09:25:41 GMT

I take the example in the link and it works on my sample dataset. One quest, how to tell whether to take 2-step regression or just 1-step regression.

When turner is 95 or 98 or 99, somewhere just 1-step (one line) is preferred ... what criteria to take to tell 2-step regression is preferred?!

%let slope1=-2; %let slope2=5; %let turner=90;

data indata;
do i=1 to 100;
	if i<&turner. then y=5*ranuni(i)+i*&slope1.;
	else y=10*ranuni(i)+i*&slope2.-&turner.*&slope2.+&turner.*&slope1.;
	output;
end;
run;quit;

ods listing gpath="%sysfunc(getoption(work))";
proc sgplot data=indata;
scatter x=i y=y; 
run;quit;

Re: How to do Regression on Broken/Turning Line?!

PaigeMiller — Fri, 16 Sep 2022 10:03:59 GMT

This is a question I have never had before, and I admit I don't know if there is a standard answer in this situation.

I assume you could compare the fits (via R-squared or RMSE or something similar) of the two models (one model with the break point and one model without). Something like the Extra Sum of Squares F-test makes sense to me. https://www.sfu.ca/~lockhart/richard/350/08_2/lectures/FTests/web.pdf

Re: How to do Regression on Broken/Turning Line?!

hellohere — Fri, 16 Sep 2022 10:23:53 GMT

My toy sample's plot and Output are below.

I bet EST/STDERR on slope1 and slope2[F Test] might be close, since F-Test between 2-step model and 1-step model is not directly available.

BUT I am note sure since not of my expertise.


                                                             The NLIN Procedure

                                                    Estimation Summary (Not Converged)

                                                   Subiterations                      98
                                                   Average Subiterations             3.5
                                                   R                            0.015325
                                                   PPC(intcpt0)                 0.019684
                                                   RPC                                 .
                                                   Object                       3.72E-16
                                                   Objective                     8.39E13
                                                   Observations Read                  66
                                                   Observations Used                  66
                                                   Observations Missing                0


                                                                   Sum of        Mean               Approx
                                 Source                    DF     Squares      Square    F Value    Pr > F

                                 Model                      3    1.124E15    3.746E14     276.85    <.0001
                                 Error                     62     8.39E13    1.353E12
                                 Corrected Total           65    1.208E15


                                                                      Approx
                                        Parameter      Estimate    Std Error    Approximate 95% Confidence Limits

                                        intcpt0         1453782       526391      401542     2506022
                                        slope1          -482416      41921.6     -566216     -398615
                                        slope2          -121843      13352.5     -148534    -95151.8
                                        x0              21.0000       1.6745     17.6526     24.3474


                                                      Approximate Correlation Matrix
                                                  intcpt0          slope1          slope2              x0

                                  intcpt0       1.0000000      -0.8760376      -0.0000000      -0.4054904
                                  slope1       -0.8760376       1.0000000       0.0000000       0.6943031
                                  slope2       -0.0000000       0.0000000       1.0000000       0.5086293
                                  x0           -0.4054904       0.6943031       0.5086293       1.0000000

Re: How to do Regression on Broken/Turning Line?!

PaigeMiller — Fri, 16 Sep 2022 10:31:28 GMT

You would have to code the equivalent of something like extra-sum-of-squares F-test yourself, based upon your PROC NLIN output. (Or do it with a calculator)

Re: How to do Regression on Broken/Turning Line?!

Ksharp — Fri, 16 Sep 2022 12:23:57 GMT

@Rick_SAS wrote a blog about this topic before.

Re: How to do Regression on Broken/Turning Line?!

Rick_SAS — Fri, 16 Sep 2022 12:43:46 GMT

This is called a piecewise linear model when each regression curve is linear. In general, these are called segmented models. See

Re: How to do Regression on Broken/Turning Line?!

PaigeMiller — Fri, 16 Sep 2022 12:48:18 GMT

Is there an index/list of topics that @Rick_SAS has written about? Or should we just assume that he has written about every possible topic 😁🙀👍 ?

Re: How to do Regression on Broken/Turning Line?!

Rick_SAS — Fri, 16 Sep 2022 13:07:55 GMT

No index, but you can search SAS blogs by using the SITE: keyword in an internet search engine. For example:

piecewise regression site:blogs.sas.com

will locate the two articles that I linked to.

Re: How to do Regression on Broken/Turning Line?!

PaigeMiller — Fri, 16 Sep 2022 13:19:48 GMT

Thanks! Actually, Google is almost as good as an index/list.

Re: How to do Regression on Broken/Turning Line?!

hellohere — Fri, 16 Sep 2022 22:22:00 GMT

The example given at R seems show that EST/STDERR for each slope is not a good tell (slope1 abs-t-value <1.0).

The visual tells a segmented line there.

My quest is how to tell whether 2-step/segment is better than 1-step (or 3-step is better than 2-step).

> slope(o)
$x
            Est.  St.Err.  t value CI(95%).l  CI(95%).u
slope1 -0.058993 0.062105 -0.94991  -0.18230 0.06431700
slope2  1.414600 0.046151 30.65100   1.32290 1.50620000
slope3 -0.142800 0.071994 -1.98350  -0.28575 0.00014734

# http://127.0.0.1:16212/library/segmented/html/segmented.html

set.seed(12)
xx<-1:100
zz<-runif(100)
yy<-2+1.5*pmax(xx-35,0)-1.5*pmax(xx-70,0)+15*pmax(zz-.5,0)+rnorm(100,0,2)
dati<-data.frame(x=xx,y=yy,z=zz)
out.lm<-lm(y~x,data=dati)

#the simplest example: the starting model includes just 1 covariate 
#.. and 1 breakpoint has to be estimated for that
o<-segmented(out.lm) #1 breakpoint for x

#the single segmented variable is not in the starting model and 1 breakpoint for that:
#... you need to specify the variable via seg.Z, but no starting value for psi
o<-segmented(out.lm, seg.Z=~z)
#note the leftmost slope is constrained to be zero (since out.lm does not include z)

#2 segmented variables, 1 breakpoint each (again no need to specify npsi or psi)
o<-segmented(out.lm,seg.Z=~z+x)


#1 segmented variable, 2 breakpoints: you have to specify starting values (vector) for psi:
o<-segmented(out.lm,seg.Z=~x,psi=c(30,60), control=seg.control(display=FALSE))

#or by specifying just the *number* of breakpoints
#o<-segmented(out.lm,seg.Z=~x, npsi=2, control=seg.control(display=FALSE)) 

slope(o) #the slopes of the segmented relationship

dat2 = data.frame(x = xx, y = broken.line(o)$fit)

ggplot(dati, aes(x = x, y = y)) +
  geom_point() +
  geom_line(data = dat2, color = 'blue')

Re: How to do Regression on Broken/Turning Line?!

hellohere — Fri, 16 Sep 2022 22:32:58 GMT

Well, EST/STDERR tells slopex=0.0 or not. But still any criteria to tell whether take n-segment rather than (n-1) segment?!

Re: How to do Regression on Broken/Turning Line?!

hellohere — Sat, 17 Sep 2022 01:11:29 GMT

I set turner=95, and fetch ESS and follow "Hypothesis Test, 3.2" and get

ESS_F=((315653-318355)/2)/(318355/96)=-0.4, somewhere must be not right ??!


%let slope1=-2; %let slope2=5; %let turner=95;

data indata;
do ind=1 to 100;
	if ind<&turner. then y=10*ranuni(ind)+ind*&slope1.;
	else y=20*ranuni(ind)+ind*&slope2.-&turner.*&slope2.+&turner.*&slope1.;
	output;
end;
run;quit;

ods listing gpath="%sysfunc(getoption(work))";
proc sgplot data=indata;
scatter x=ind y=y; 
run;quit;

%let inds=indata;
%let xvar=ind; %let yvar=y;

	proc nlin data=&inds.;
	   parms intcpt0=0 slope1=0 slope2=0 x0=&st_x0;
	   if (ind < x0) then
	        mean = intcpt0+ind*slope1;
	   else mean = ind*slope2-x0*slope2+x0*slope1+intcpt0;
	   model &yvar. = mean;

	   output out=&inds._sel_out predicted=&yvar._p;
	run;

	proc reg data=&inds.;
	model y=ind;run;quit;


                                                             The REG Procedure
                                                               Model: MODEL1
                                                           Dependent Variable: y

                                                  Number of Observations Read         100
                                                  Number of Observations Used         100


                                                            Analysis of Variance

                                                                   Sum of           Mean
                               Source                   DF        Squares         Square    F Value    Pr > F

                               Model                     1         315653         315653    8832.21    <.0001
                               Error                    98     3502.40844       35.73886
                               Corrected Total          99         319156


                                            Root MSE              5.97820    R-Square     0.9890
                                            Dependent Mean      -94.98238    Adj R-Sq     0.9889
                                            Coeff Var            -6.29401


                                                            Parameter Estimates

                                                         Parameter       Standard
                                    Variable     DF       Estimate          Error    t Value    Pr > |t|

                                    Intercept     1        3.30749        1.20466       2.75      0.0072
                                    ind           1       -1.94633        0.02071     -93.98      <.0001


                                                             The NLIN Procedure

                                                                   Sum of        Mean               Approx
                                 Source                    DF     Squares      Square    F Value    Pr > F

                                 Model                      3      318355      106118    12731.3    <.0001
                                 Error                     96       800.2      8.3353
                                 Corrected Total           99      319156


                                                                      Approx
                                        Parameter      Estimate    Std Error    Approximate 95% Confidence Limits

                                        intcpt0          5.6447       0.6003      4.4530      6.8363
                                        slope1          -2.0169       0.0110     -2.0387     -1.9951
                                        slope2           4.5742       0.6901      3.2043      5.9441
                                        x0              94.3985       0.3816     93.6411     95.1559


                                                      Approximate Correlation Matrix
                                                  intcpt0          slope1          slope2              x0

                                  intcpt0       1.0000000      -0.8683135      -0.0000000      -0.1189745
                                  slope1       -0.8683135       1.0000000       0.0000000       0.2046574
                                  slope2       -0.0000000       0.0000000       1.0000000       0.8511409
                                  x0           -0.1189745       0.2046574       0.8511409       1.0000000

Re: How to do Regression on Broken/Turning Line?!

Rick_SAS — Sat, 17 Sep 2022 10:13:04 GMT

The usual way to compare two different regression models is to compare some criteria such as the AIC, BIC, or SBC. These criteria are not produced automatically by PROC NLIN, but you can use PROC NLMIXED (which has a similar syntax) to get them. Look at the "Fit Statistics" table.

Re: How to do Regression on Broken/Turning Line?!

hellohere — Sun, 18 Sep 2022 00:58:42 GMT

Would you give out an example with dataset or a link?! Dear Paige gives out a great link on hypo test between

models, https://www.sfu.ca/~lockhart/richard/350/08_2/lectures/FTests/web.pdf

Somehow, when I run my toy dataset, the F-test shows up confusing.

%let slope1=-2; %let slope2=5; %let turner=95;

data indata;
do ind=1 to 100;
	if ind<&turner. then y=10*ranuni(ind)+ind*&slope1.;
	else y=20*ranuni(ind)+ind*&slope2.-&turner.*&slope2.+&turner.*&slope1.;
	output;
end;
run;quit;

ods listing gpath="%sysfunc(getoption(work))";
proc sgplot data=indata;
scatter x=ind y=y; 
run;quit;

%let inds=indata;
%let xvar=ind; %let yvar=y;

	proc nlin data=&inds.;
	   parms intcpt0=0 slope1=0 slope2=0 x0=&st_x0;
	   if (ind < x0) then
	        mean = intcpt0+ind*slope1;
	   else mean = ind*slope2-x0*slope2+x0*slope1+intcpt0;
	   model &yvar. = mean;

	   output out=&inds._sel_out predicted=&yvar._p;
	run;

	proc reg data=&inds.;
	model y=ind;run;quit;