Nonlinear Models Action Set

Segmented Model

This example illustrates the fitting of a segmented model by using a hypothetical data set. Suppose you are interested in fitting a model that consists of two segments that connect in a smooth fashion. For example, the following model states that the mean of Y is a quadratic function in x for values of x less than and that the mean of Y is constant for values of x greater than :

normal upper E left-bracket upper Y vertical-bar x right-bracket equals StartLayout Enlarged left-brace 1st Row 1st Column alpha plus beta x plus gamma x squared 2nd Column normal i normal f x less-than x 0 2nd Row 1st Column c 2nd Column normal i normal f x greater-than-or-equal-to x 0 EndLayout

In this model equation, , , and are the coefficients of the quadratic segment, and c is the plateau of the mean function. The nlmod action can fit such a segmented model even when the join point, , is unknown.

Suppose you also want to impose conditions on the two segments of the model. First, the curve should be continuous—that is, the quadratic segment and the plateau section need to meet at . Second, the curve should be smooth—that is, the first derivatives of the two segments with respect to x need to coincide at .

The continuity condition requires that

c equals normal upper E left-bracket upper Y vertical-bar x 0 right-bracket equals alpha plus beta x 0 plus gamma x 0 squared

The smoothness condition requires that

StartFraction partial-differential normal upper E left-bracket upper Y vertical-bar x 0 right-bracket Over partial-differential x EndFraction equals beta plus 2 gamma x 0 identical-to 0

If you solve for and substitute your result in the expression for c, the two conditions jointly imply that

StartLayout 1st Row 1st Column x 0 2nd Column equals negative beta slash 2 gamma 2nd Row 1st Column c 2nd Column equals alpha minus beta squared slash 4 gamma EndLayout

Although there are five unknowns, the model contains only three independent parameters. Together the continuity and smoothness restrictions completely determine two parameters, given the other three.

The following DATA step creates the data table for this example in your CAS session:

data mycas.a;
   input y x @@;
   datalines;
.46 1  .47  2 .57  3 .61  4 .62  5 .68  6 .69  7
.78 8  .70  9 .74 10 .77 11 .78 12 .74 13 .80 13
.80 15 .78 16
 ;

The following nlmod action program fits the segmented model to the data table a:

proc cas;
   nonlinear.nlmod
      table='a',
      parms={{name='alpha',vals=.45},
             {name='beta',vals=.05},
             {name='gamma',vals=-.0025}},
      nlmodCode='x0 = -.5*beta / gamma;
                 if (x < x0) then
                   yp = alpha + beta*x  + gamma*x*x;
                 else
                   yp = alpha + beta*x0 + gamma*x0*x0;',
      model={depvar='y',
             dist='residual',
             distparms='yp'},
      estimate={{label='join point', expression='-beta/2/gamma'},
                {label='plateau value c', expression='alpha-beta**2/(4*gamma)'}},
      predout={name='b'},
      predict={{label='predicted', expression='yp', pred='yp'}},
      id={'y','x'};
run;

The table parameter names the input data table. The parameters of the model are , , and , which are represented in the nlmod action parameters by the variables alpha, beta, and gamma, respectively. The parms parameter specifies the initial values for these model parameters. To model the two segments, a conditional statement assigns the appropriate expression to the mean function, depending on the value of . You use the nlmodCode parameter to specify this conditional assignment. The estimate parameter computes the values of and c. The predict parameter computes predicted values for plotting and saves them to data table b, which is specified in the name subparameter of the predout parameter. Along with the predicted values, the data set variables y and x are also saved to data table b; they are specified using the id parameter.

The results from fitting this model are shown in Figure 10 through Figure 12. The iterative optimization converges after six iterations (Figure 10). Figure 11 shows the estimated parameters. Figure 12 indicates that the join point is 12.7477 and the plateau value is 0.7775.

Figure 10: Nonlinear Least Squares Iterative Phase

Results from nonlinear.nlmod

Iteration History
Iteration	Evaluations	Objective Function	Change	Maximum Gradient
0	5	0.0035144531		7.184063
1	2	0.0007352716	0.00277918	2.145337
2	2	0.0006292751	0.00010600	0.032551
3	2	0.0006291261	0.00000015	0.002952
4	2	0.0006291244	0.00000000	0.000238
5	2	0.0006291244	0.00000000	0.000023
6	2	0.0006291244	0.00000000	2.313E-6

Convergence criterion (GCONV=1E-8) satisfied.

Figure 11: Least Squares Analysis of the Quadratic Model

Analysis of Variance
Source	DF	Sum of Squares	Mean Square	F Value	Approx Pr > F
Model	2	0.1769	0.0884	114.22	<.0001
Error	13	0.0101	0.000774
Corrected Total	15	0.1869

Parameter Estimates
Parameter	Estimate	Standard Error	DF	t Value	Approx Pr > \|t\|	Approximate 95% Confidence Limits
Parameter	Estimate	Standard Error	DF	t Value	Approx Pr > \|t\|	Lower	Upper
alpha	0.3921	0.0267	13	14.70	<.0001	0.3345	0.4497
beta	0.0605	0.00842	13	7.18	<.0001	0.0423	0.0787
gamma	-0.00237	0.000551	13	-4.30	0.0009	-0.00356	-0.00118

Figure 12: Additional Estimates for the Quadratic Model

Additional Estimates
Label	Estimate	Standard Error	DF	t Value	Approx Pr > \|t\|	Alpha	Approximate Confidence Limits
Label	Estimate	Standard Error	DF	t Value	Approx Pr > \|t\|	Alpha	Lower	Upper
join point	12.7477	1.2781	13	9.97	<.0001	0.05	9.9864	15.5089
plateau value c	0.7775	0.0123	13	63.11	<.0001	0.05	0.7509	0.8041

By using the predicted values and the y values from data table b, you can plot the observed versus predicted values to check the goodness of fit of the model.

Segmented Model

This section contains Lua code for the analysis in the CASL version of this example, which contains details about the results.

Note: In order to run this code, the data that are described in the CASL version need to be accessible to the CAS server. One way to do this is to convert the a data to the comma-separated-value (CSV) file a.csv and then use the following code to load the CSV file into CAS:

s:loadtable{casLib="casuser", path="a.csv"}

For more information about coding in Lua, see Getting Started with SAS Viya for Lua and SAS Viya: System Programming Guide.

The following code loads the nonlinear action set and then uses the nlmod action in the nonlinear action set to fit a segmented model to the data table a:

s:loadactionset{actionset='nonlinear'}
m=s:nlmod{table={name='a'},
      parms={{name='alpha',vals=.45},
             {name='beta',vals=.05},
             {name='gamma',vals=-.0025}},
      nlmodCode='x0 = -.5*beta / gamma; \z
                 if (x < x0) then \z
                   yp = alpha + beta*x  + gamma*x*x; \z
                 else \z
                   yp = alpha + beta*x0 + gamma*x0*x0;',
      model={depvar='y',
             dist='residual',
             distparms='yp'},
      estimate={{label='join point', expression='-beta/2/gamma'},
                {label='plateau value c', expression='alpha-beta**2/(4*gamma)'}},
      predout={name='b'},
      predict={{label='predicted', expression='yp', pred='yp'}},
      id={'y','x'} }

The following commands display the tables that are produced by this action call:

print(m.IterHistory)
print(m.ConvergenceStatus)
print(m.ANOVA)
print(m.ParameterEstimates)
print(m.AdditionalEstimates)
print(m.Timing)

For details about the results of this analysis, see the CASL version of this example.

Segmented Model

This section contains Python code for the analysis in the CASL version of this example, which contains details about the results.

s.upload_file('a.csv')

For more information about coding in Python, see Getting Started with SAS Viya for Python and SAS Viya: System Programming Guide.

The following code loads the nonlinear action set and then uses the nlmod action in the nonlinear action set to fit a segmented model to the data table a:

s.loadactionset(actionset='nonlinear')
m=s.nlmod(table='a',
      parms=[{'name':'alpha','vals':.45},
             {'name':'beta','vals':.05},
             {'name':'gamma','vals':-.0025}],
      nlmodCode="""x0 = -.5*beta / gamma;
                 if (x < x0) then
                   yp = alpha + beta*x  + gamma*x*x;
                 else
                   yp = alpha + beta*x0 + gamma*x0*x0;""",
      model={'depvar':'y',
             'dist':'residual',
             'distparms':'yp'},
      estimate=[{'label':'join point', 'expression':'-beta/2/gamma'},
                {'label':'plateau value c', 'expression':'alpha-beta**2/(4*gamma)'}],
      predout={'name':'b'},
      predict=[{'label':'predicted', 'expression':'yp', 'pred':'yp'}],
      id=['y','x'] )

The following commands display the tables that are produced by this action call:

print(m.IterHistory)
print(m.ConvergenceStatus)
print(m.ANOVA)
print(m.ParameterEstimates)
print(m.AdditionalEstimates)
print(m.Timing)

The following command displays all the tables that are produced by this action call along with their full pathnames:

print(m)

For details about the results of this analysis, see the CASL version of this example.

Segmented Model

This section contains R code for the analysis in the CASL version of this example, which contains details about the results.

m <- cas.read.csv(s, "a.csv", casOut=list(name="a"))

For more information about coding in R, see Getting Started with SAS Viya for R and SAS Viya: System Programming Guide.

The following code loads the nonlinear action set and then uses the nlmod action in the nonlinear action set to fit a segmented model to the data table a:

m <- cas.builtins.loadActionSet(s, actionset='nonlinear')
m <- cas.nonlinear.nlmod(s,
      table='a',
      parms=list(list(name='alpha',vals=.45),
             list(name='beta',vals=.05),
             list(name='gamma',vals=-.0025)),
      nlmodCode='x0 = -.5*beta / gamma;
                 if (x < x0) then
                   yp = alpha + beta*x  + gamma*x*x;
                 else
                   yp = alpha + beta*x0 + gamma*x0*x0;',
      model=list(depvar='y',
             dist='residual',
             distparms='yp'),
      estimate=list(list(label='join point', expression='-beta/2/gamma'),
                list(label='plateau value c', expression='alpha-beta**2/(4*gamma)')),
      predout=list(name='b'),
      predict=list(list(label='predicted', expression='yp', pred='yp')),
      id=list('y','x') )

The following commands display the tables that are produced by this action call:

print(m$IterHistory)
print(m$ConvergenceStatus)
print(m$ANOVA)
print(m$ParameterEstimates)
print(m$AdditionalEstimates)
print(m$Timing)

The following command displays all the tables that are produced by this action call along with their full pathnames:

print(m)

For details about the results of this analysis, see the CASL version of this example.

Last updated: September 13, 2022