topic Re: Estimate inflection point for 2 linear lines in SAS Programming

Estimate inflection point for 2 linear lines

PamG — Mon, 21 Aug 2023 19:01:18 GMT

I have a data where the slope changes at a certain point. How do you find that point of inflection for the two linear lines? I think PROC NLIN would work but I can't get it to run.

Sample code is given below as I cannot share my data. The point of inflection is 50.

DATA data1;
do i=1 to 1000;
x=rand("Integer", 1, 50);
Y=3 + 6*x + rannor(345);
drop i;
OUTPUT;
end;
RUN;

DATA data2;
do i=1 to 1000;
x=rand("Integer", 50, 100);
Y=2 + 8*x + rannor(345);
drop i;
OUTPUT;
end;
RUN;

data full; set data1 data2;run;

proc nlin data=full;
   parms alpha1=.45 beta1=.05 alpha2=.45 beta2=.05;
x0=5*beta1/alpha1;
 model y = (x <x0)*(alpha1+beta1*x) +
             (x>=x0)*(alpha2+beta2*x);
FILE log;put  / x0=  ; 
output out=b predicted=yp;
run;

Re: Estimate inflection point for 2 linear lines

PaigeMiller — Mon, 21 Aug 2023 19:05:02 GMT

I think PROC NLIN would work but I can't get it to run.

What is wrong with your code?

Re: Estimate inflection point for 2 linear lines

PamG — Mon, 21 Aug 2023 19:15:18 GMT

The Hessian matrix is singular so I do not get parameter estimates. I am not sure the code I have written for linear model is correct.

Re: Estimate inflection point for 2 linear lines

PaigeMiller — Mon, 21 Aug 2023 19:39:29 GMT

@PamG wrote:

The Hessian matrix is singular so I do not get parameter estimates. I am not sure the code I have written for linear model is correct.

Thank you. In the future, key information like this should be in your first message.

This discusses the error message that says the Hessian is singular: https://communities.sas.com/t5/Statistical-Procedures/Problem-in-Proc-NLIN-The-approximate-Hessian-is-singular/td-p/136949 There may be other similar discussions here in the SAS communities, and you should search for them.

Re: Estimate inflection point for 2 linear lines

sbxkoenk — Mon, 21 Aug 2023 19:49:37 GMT

Hello @PamG ,

If your only interest is to know about the inflection point ,

then you can use PROC UCM (SAS/ETS) for structural break detection in a univariate time series.

I use your variable X as a time indicator. This is possible because in your sample data the X-values are equally spaced (fixed width interval).
Not sure how close X is to a real time series in your actual data (?).

DATA data1;
do i=1 to 1000;
x=rand("Integer", 1, 50);
Y=3 + 6*x + rannor(345);
drop i;
OUTPUT;
end;
RUN;

DATA data2;
do i=1 to 1000;
x=rand("Integer", 50, 100);
Y=2 + 8*x + rannor(345);
drop i;
OUTPUT;
end;
RUN;

data full; set data1 data2; run;
proc sort data=full; by x; run;

proc sgplot data=full;
  scatter x=x y=y;
run;

/* Accumulating Transactional Data into Time Series Data   */
/* Goal is to get Only 1 observation per unique value of X */
proc timeseries data=full out=tsfull;
   id x         interval=day
                accumulate=median
                setmiss=0
                ;
   var y;
run;

proc sgplot data=tsfull;
  scatter x=x y=y;
run;

ods trace off;
/* Detection of Level Shift by unobserved components model (UCM).                               */
/* (The UCMs are also called structural models in the time series literature.)                  */
/* You can check for breaks in the level by using the CHECKBREAK option in the LEVEL statement. */ 
/* The following statements fit a simple locally constant level plus constant slope             */
/* plus error model to the series:                                                              */ 
ODS EXCLUDE EstimationSpan ForecastSpan InitialParameters FitSummary
            ParameterEstimates FitStatistics ComponentSignificance
            TrendInformation Forecasts ForecastsOnlyPlot ;
ODS OUTPUT OutlierSummary=work.OutlierSummary;
proc ucm data=tsfull;
   id x interval=day;
   model y;
   irregular;
/* The fitted trend is nearly linear,                                          */ 
/* So I fix fix the level and/or slope variance components to zero             */
/* (via the VARIANCE=0 and NOEST options on the LEVEL and/or SLOPE statements) */
   level variance=0 noest /* plot=smooth */ checkbreak;
   slope variance=0 noest;
  *estimate;
  *forecast plot=decomp;
run;

data _NULL_;
 set work.OutlierSummary;
 file print;
 put 'inflection point is at ' time= 12.;
run;
/* end of program */

BR,

Koen

Re: Estimate inflection point for 2 linear lines

sbxkoenk — Mon, 21 Aug 2023 19:51:33 GMT

On top of PROC UCM post above ,

Re: Estimate inflection point for 2 linear lines

ballardw — Mon, 21 Aug 2023 20:24:03 GMT

@PamG wrote:

I have a data where the slope changes at a certain point. How do you find that point of inflection for the two linear lines?

"Linear Lines" have no inflection point where the slope changes. Curves may.

Are you assuming both of your curves have inflection points at the same value? That may be pretty broad assumption.

Are you using a different definition of inflection other than a point on the curve where the slope changes sign?

BTW, if your data is so sensitive that a group of x,y pairs with no other identification causes concerns then perhaps you shouldn't be asking any questions on an open forum.

Re: Estimate inflection point for 2 linear lines

PamG — Mon, 21 Aug 2023 20:50:39 GMT

Thanks Koen. I do not have a time series data and I do not have a jump in the line. I am so sorry but I do not know how to create sample datasets to show continuous lines but with different slopes. I have (y,x) points and I would like to know where the joint point is that has no breaks.

Re: Estimate inflection point for 2 linear lines

Ksharp — Tue, 22 Aug 2023 11:37:17 GMT

https://blogs.sas.com/content/iml/2020/12/14/segmented-regression-sas.html
https://blogs.sas.com/content/iml/2017/04/05/nonsmooth-models-spline-effects.html

Re: Estimate inflection point for 2 linear lines

FreelanceReinh — Tue, 22 Aug 2023 18:38:20 GMT

Hello @PamG,

Another option: Try PROC FASTCLUS to get a first estimate of the point of inflection and then use PROC ROBUSTREG for fine-tuning (see the proportions of outliers in the "Diagnostic Summary" table of the output) and estimation of the regression coefficients.

proc fastclus data=full out=clust maxclusters=2 maxiter=100;
var x y;
run;

proc means data=clust;
class cluster;
var x;
run;

proc sort data=clust;
by cluster;
run;

proc robustreg data=clust method=lts seed=2718;
by cluster;
model y=x;
run;