Interaction term in non-linear (nl) regression models

[Peta Hitchens]

Good afternoon All

Sample data is supplied at the bottom of this post.

I have an equation defined below that I have modelled using nl (nonlinear least-squares estimation).
y = {b0}*sqrt(1-x)*({b1}*x+1-{b1}

nl (y= {b0}*sqrt(1-x)*({b1}*x+1-{b1})), nolog variables(x)
margins, at(x=(0.95(.0025)1))
marginsplot
addplot: scatter y x, jitter(1) msym(oh) mcolor(gs10)

I have been able to graph the differences between binary (e.g. sex) and continuous (e.g. age) predictors , but I would like to be able to determine whether there is a significant difference between the two curves (this would be easy if it was a quadratic fit for example: regress y c.x##c.x##i.sex).

*determine equaton fit by sex, combine curves
set more off
nl (y= {b0}*sqrt(1-x)*({b1}*x+1-{b1})) if sex==0, nolog variables(x)
predict int_female
nl (y= {b0}*sqrt(1-x)*({b1}*x+1-{b1})) if sex==1, nolog variables(x)
predict int_male
scatter y x || lowess int_female x || lowess int_male x, sort

My attempt at including interaction terms (main effect for x and sex, and interaction term between them) in the equation is below, but I am not at all certain whether this is correct at all! The p-value from this sample data suggests that there is no significant difference between sexes (p=0.619):

nl (y = {b0}*sqrt(1-x)*({b1}*x+1-{b1})+{b2}*{sex}+{i1}*({b0}*sqrt(1-x)*({b1}*x+1-{b1}))*sex), nolog variables(x sex)
margins, over(sex) at(x=(0.95(.005)1))
marginsplot

I'd be super grateful if someone could confirm whether or not the interaction term is correctly fitted in the nl equation!

Thanks for any help or advice,
/Peta

------------------------------------------------------------------------------------------------------------------------------------

* Example generated by -dataex-. To install: ssc install dataex
clear
input float y double x byte(sex age)
1.52796 .9892868 1 2
1.89666 .9797394 0 6
1.96257 .978637 1 3
1.98382 .9780292 0 4
2.04352 .9726008 1 3
2.09717 .9673443 1 4
2.1616 .9697651 1 2
2.35337 .9626659 0 4
2.36128 .9639572 0 3
2.39121 .9704398000000001 1 7
2.39734 .9699420999999999 1 6
2.42883 .9548707000000001 1 7
2.4722 .9531984 1 2
2.47434 .9673583000000001 1 3
2.48315 .9524335 1 2
2.54352 .952538 0 3
2.55237 .9584881 0 3
2.62491 .9660379 1 7
2.65316 .9530894000000001 1 4
2.66142 .9520308 1 4
end
label values sex sex
label def sex 0 "female", modify
label def sex 1 "male", modify

[Joseph Coveney]

What you show doesn't really converge to anything interesting when I try it on your sample dataset (see below, first model fitted).

You might want to consider something along these lines:

nl ( y = {b0} * sqrt(1-x) * ( {b1} * x + 1 - {b1} ) + ///
    {s} * sex + ///
    {sb0} * sqrt(1-x*sex) * ( {sb1} * x * sex + 1 - {sb1} ) )

which is illustrated below. Because it has five parameter estimates with only 20 observations, I first fitted the female-only (sex == 0) subset of data in order to get some starting values, which helps stability and speeds up convergence.

.ÿversionÿ14.2

.ÿ
.ÿclearÿ*

.ÿsetÿmoreÿoff

.ÿ
.ÿinputÿfloatÿyÿdoubleÿxÿbyte(sexÿage)

ÿÿÿÿÿÿÿÿÿÿÿÿÿyÿÿÿÿÿÿÿÿÿÿÿxÿÿÿÿÿÿÿsexÿÿÿÿÿÿÿage
ÿÿ1.ÿ1.52796ÿ.9892868ÿ1ÿ2
ÿÿ2.ÿ1.89666ÿ.9797394ÿ0ÿ6
ÿÿ3.ÿ1.96257ÿ.978637ÿ1ÿ3
ÿÿ4.ÿ1.98382ÿ.9780292ÿ0ÿ4
ÿÿ5.ÿ2.04352ÿ.9726008ÿ1ÿ3
ÿÿ6.ÿ2.09717ÿ.9673443ÿ1ÿ4
ÿÿ7.ÿ2.1616ÿ.9697651ÿ1ÿ2
ÿÿ8.ÿ2.35337ÿ.9626659ÿ0ÿ4
ÿÿ9.ÿ2.36128ÿ.9639572ÿ0ÿ3
ÿ10.ÿ2.39121ÿ.9704398000000001ÿ1ÿ7
ÿ11.ÿ2.39734ÿ.9699420999999999ÿ1ÿ6
ÿ12.ÿ2.42883ÿ.9548707000000001ÿ1ÿ7
ÿ13.ÿ2.4722ÿ.9531984ÿ1ÿ2
ÿ14.ÿ2.47434ÿ.9673583000000001ÿ1ÿ3
ÿ15.ÿ2.48315ÿ.9524335ÿ1ÿ2
ÿ16.ÿ2.54352ÿ.952538ÿ0ÿ3
ÿ17.ÿ2.55237ÿ.9584881ÿ0ÿ3
ÿ18.ÿ2.62491ÿ.9660379ÿ1ÿ7
ÿ19.ÿ2.65316ÿ.9530894000000001ÿ1ÿ4
ÿ20.ÿ2.66142ÿ.9520308ÿ1ÿ4
ÿ21.ÿend

.ÿ
.ÿ*
.ÿ*ÿVerbatimÿfromÿyourÿpost
.ÿ*
.ÿnlÿ(yÿ=ÿ{b0}*sqrt(1-x)*({b1}*x+1-{b1})+{b2}*{sex}+{i1}*({b0}*sqrt(1-x)*({b1}*x+1-{b1}))*sex),ÿnologÿvariables(xÿsex)
(obsÿ=ÿ20)


ÿÿÿÿÿÿSourceÿ|ÿÿÿÿÿÿSSÿÿÿÿÿÿÿÿÿÿÿÿdfÿÿÿÿÿÿÿMS
-------------+----------------------------------ÿÿÿÿNumberÿofÿobsÿ=ÿÿÿÿÿÿÿÿÿ20
ÿÿÿÿÿÿÿModelÿ|ÿÿ1.4403774ÿÿÿÿÿÿÿÿÿÿ2ÿÿ.720188711ÿÿÿÿR-squaredÿÿÿÿÿ=ÿÿÿÿÿ0.8390
ÿÿÿÿResidualÿ|ÿÿ.27650116ÿÿÿÿÿÿÿÿÿ17ÿÿ.016264774ÿÿÿÿAdjÿR-squaredÿ=ÿÿÿÿÿ0.8200
-------------+----------------------------------ÿÿÿÿRootÿMSEÿÿÿÿÿÿ=ÿÿÿ.1275334
ÿÿÿÿÿÿÿTotalÿ|ÿÿ1.7168786ÿÿÿÿÿÿÿÿÿ19ÿÿ.090362031ÿÿÿÿRes.ÿdev.ÿÿÿÿÿ=ÿÿ-28.86791

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿyÿ|ÿÿÿÿÿÿCoef.ÿÿÿStd.ÿErr.ÿÿÿÿÿÿtÿÿÿÿP>|t|ÿÿÿÿÿ[95%ÿConf.ÿInterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿ/b0ÿ|ÿÿÿ15.08924ÿÿÿ.6806321ÿÿÿÿ22.17ÿÿÿ0.000ÿÿÿÿÿ13.65323ÿÿÿÿ16.52525
ÿÿÿÿÿÿÿÿÿ/b1ÿ|ÿÿÿÿ4.79515ÿÿÿ.8816482ÿÿÿÿÿ5.44ÿÿÿ0.000ÿÿÿÿÿ2.935035ÿÿÿÿ6.655266
ÿÿÿÿÿÿÿÿÿ/b2ÿ|ÿÿÿÿÿÿÿÿÿÿ0ÿÿ(constrained)
ÿÿÿÿÿÿÿÿ/sexÿ|ÿÿÿÿÿÿÿÿÿÿ0ÿÿ(constrained)
ÿÿÿÿÿÿÿÿÿ/i1ÿ|ÿÿÿÿÿ.01386ÿÿÿÿ.027389ÿÿÿÿÿ0.51ÿÿÿ0.619ÿÿÿÿ-.0439256ÿÿÿÿ.0716457
------------------------------------------------------------------------------
ÿÿParameterÿsexÿtakenÿasÿconstantÿtermÿinÿmodelÿ&ÿANOVAÿtable

.ÿ
.ÿ*
.ÿ*ÿSuggested
.ÿ*
.ÿquietlyÿnlÿ(y=ÿ{b0}*sqrt(1-x)*({b1}*x+1-{b1}))ÿifÿsex==0

.ÿlocalÿb0ÿ=ÿ_b[/b0]

.ÿlocalÿb1ÿ=ÿ_b[/b1]

.ÿ
.ÿnlÿ(ÿyÿ=ÿ{b0}ÿ*ÿsqrt(1-x)ÿ*ÿ(ÿ{b1}ÿ*ÿxÿ+ÿ1ÿ-ÿ{b1})ÿ+ÿ///
>ÿÿÿÿÿÿÿÿÿ{s}ÿ*ÿsexÿ+ÿ///
>ÿÿÿÿÿÿÿÿÿ{sb0}ÿ*ÿsqrt(1-x*sex)ÿ*ÿ(ÿ{sb1}ÿ*ÿxÿ*ÿsexÿ+ÿ1ÿ-ÿ{sb1})ÿ),ÿ///
>ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿinit(b0ÿ`b0'ÿb1ÿ`b1')ÿvariables(xÿsex)ÿnolog
(obsÿ=ÿ20)


ÿÿÿÿÿÿSourceÿ|ÿÿÿÿÿÿSSÿÿÿÿÿÿÿÿÿÿÿÿdfÿÿÿÿÿÿÿMS
-------------+----------------------------------ÿÿÿÿNumberÿofÿobsÿ=ÿÿÿÿÿÿÿÿÿ20
ÿÿÿÿÿÿÿModelÿ|ÿÿ107.56911ÿÿÿÿÿÿÿÿÿÿ5ÿÿ21.5138229ÿÿÿÿR-squaredÿÿÿÿÿ=ÿÿÿÿÿ0.9975
ÿÿÿÿResidualÿ|ÿÿ.27185283ÿÿÿÿÿÿÿÿÿ15ÿÿ.018123522ÿÿÿÿAdjÿR-squaredÿ=ÿÿÿÿÿ0.9966
-------------+----------------------------------ÿÿÿÿRootÿMSEÿÿÿÿÿÿ=ÿÿÿ.1346236
ÿÿÿÿÿÿÿTotalÿ|ÿÿ107.84097ÿÿÿÿÿÿÿÿÿ20ÿÿ5.39204837ÿÿÿÿRes.ÿdev.ÿÿÿÿÿ=ÿÿ-29.20699

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿÿÿyÿ|ÿÿÿÿÿÿCoef.ÿÿÿStd.ÿErr.ÿÿÿÿÿÿtÿÿÿÿP>|t|ÿÿÿÿÿ[95%ÿConf.ÿInterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿ/b0ÿ|ÿÿÿ17.14785ÿÿÿ5.441886ÿÿÿÿÿ3.15ÿÿÿ0.007ÿÿÿÿÿ5.548742ÿÿÿÿ28.74695
ÿÿÿÿÿÿÿÿÿ/b1ÿ|ÿÿÿ4.972012ÿÿÿ1.694938ÿÿÿÿÿ2.93ÿÿÿ0.010ÿÿÿÿÿ1.359336ÿÿÿÿ8.584688
ÿÿÿÿÿÿÿÿÿÿ/sÿ|ÿÿ-.0693768ÿÿÿ.5547121ÿÿÿÿ-0.13ÿÿÿ0.902ÿÿÿÿ-1.251718ÿÿÿÿ1.112964
ÿÿÿÿÿÿÿÿ/sb0ÿ|ÿÿ-1.062492ÿÿÿ2.499141ÿÿÿÿ-0.43ÿÿÿ0.677ÿÿÿÿ-6.389286ÿÿÿÿ4.264301
ÿÿÿÿÿÿÿÿ/sb1ÿ|ÿÿÿ.7230896ÿÿÿ.5987996ÿÿÿÿÿ1.21ÿÿÿ0.246ÿÿÿÿ-.5532215ÿÿÿÿ1.999401
------------------------------------------------------------------------------

.ÿ
.ÿtestÿ_b[/sb0]ÿ=ÿ0,ÿnotest

ÿ(ÿ1)ÿÿ[sb0]_consÿ=ÿ0

.ÿtestÿ_b[/sb1]ÿ=ÿ0,ÿaccumulate

ÿ(ÿ1)ÿÿ[sb0]_consÿ=ÿ0
ÿ(ÿ2)ÿÿ[sb1]_consÿ=ÿ0

ÿÿÿÿÿÿÿF(ÿÿ2,ÿÿÿÿ15)ÿ=ÿÿÿÿ0.74
ÿÿÿÿÿÿÿÿÿÿÿÿProbÿ>ÿFÿ=ÿÿÿÿ0.4933

.ÿ
.ÿquietlyÿmargins,ÿover(sex)ÿat(x=(0.95(.005)1))

.ÿmarginsplotÿ,ÿxlabel(0.95(0.01)1)ÿylabel(ÿ,ÿangle(horizontal)ÿnogrid)ÿ///
>ÿÿÿÿÿÿÿÿÿplot1opts(msymbol(none)ÿlcolor(red))ÿplot2opts(msymbol(none)ÿlcolor(blue))ÿ///
>ÿÿÿÿÿÿÿÿÿnociÿtitle("")ÿlegend(off)

ÿÿVariablesÿthatÿuniquelyÿidentifyÿmargins:ÿxÿsex

.ÿ
.ÿexit

endÿofÿdo-file


.

[Peta Hitchens]

Thanks Joseph.
I have tried your suggestion on my full dataset (n=210; full dataset at end of post), but surprised that it suggests a significant difference between males and females (p<0.001), if I am interpreting that correctly?

Though coefficients sb0 and sb1 are both non-sig at p>0.5. I would not expect significance by looking at the graphs and because the results of a quadratic equation with sex interaction is not sig either.

quietly nl (y= {b0}*sqrt(1-x)*({b1}*x+1-{b1})) if sex==0
local b0 = _b[/b0]
local b1 = _b[/b1]
nl (y = {b0}*sqrt(1-x)*({b1}*x+1-{b1}) + {s}*sex + {sb0}*sqrt(1-x*sex)*({sb1}*x*sex+1-{sb1})), initial(b0 `b0' b1 `b1') variables(x sex) nolog

Results:

      Source |      SS            df       MS
-------------+----------------------------------    Number of obs =        210
       Model |  2315.8545          5  463.170903    R-squared     =     0.9909
    Residual |  21.209932        205  .103463083    Adj R-squared =     0.9907
-------------+----------------------------------    Root MSE      =   .3216568
       Total |  2337.0644        210  11.1288783    Res. dev.     =   114.5002

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         /b0 |   8.665484   .6863062    12.63   0.000     7.312361    10.01861
         /b1 |    .120811   .1864269     0.65   0.518    -.2467489     .488371
          /s |   .6708679    .128342     5.23   0.000     .4178283    .9239074
        /sb0 |   .3284166   .6028605     0.54   0.587    -.8601853    1.517019
        /sb1 |  -1.296689   3.936345    -0.33   0.742    -9.057601    6.464223
------------------------------------------------------------------------------

quietly margins, over(sex) at(x=(0.65(0.05)1))
marginsplot
test _b[/sb0] = 0
test _b[/sb1] = 0, accumulate

Results of significance test of interaction (I think):

test _b[/sb0] = 0
 ( 1)  [sb0]_cons = 0

       F(  1,   205) =    0.30
            Prob > F =    0.5865

test _b[/sb1] = 0, accumulate

 ( 1)  [sb0]_cons = 0
 ( 2)  [sb1]_cons = 0

       F(  2,   205) =   40.38
            Prob > F =    0.0000

Thanks
/Peta

* Example generated by -dataex-. To install: ssc install dataex
clear
input float y double x byte(sex age)
  1.52796          .9892868 1 2
  1.89666          .9797394 0 6
  1.96257           .978637 1 3
  1.98382          .9780292 0 4
  2.04352          .9726008 1 3
  2.09717          .9673443 1 4
   2.1616          .9697651 1 2
  2.35337          .9626659 0 4
  2.36128          .9639572 0 3
  2.39121 .9704398000000001 1 7
  2.39734 .9699420999999999 1 6
  2.42883 .9548707000000001 1 7
   2.4722          .9531984 1 2
  2.47434 .9673583000000001 1 3
  2.48315          .9524335 1 2
  2.54352           .952538 0 3
  2.55237          .9584881 0 3
  2.60276          .8288718 1 3
  2.62491          .9660379 1 7
  2.65316 .9530894000000001 1 4
  2.66142          .9520308 1 4
  2.77439          .9427006 1 5
  2.82644 .9328816000000001 1 3
  2.86415          .9141257 1 3
  2.88059 .8943971999999999 0 6
  2.88923 .9511056999999999 1 2
   2.9161 .9327715000000001 1 4
  2.94116 .9337726000000001 1 3
  2.95274          .9241852 1 3
  2.95658          .9324329 1 6
  2.97797          .9354988 0 2
  3.05192 .9464950000000001 1 2
 1.394259  .990999984741211 1 3
  3.07709           .932842 1 4
  3.09643          .9204298 0 2
  3.12332 .9248816999999999 1 4
   3.1273 .9162184000000001 0 2
  3.15498          .9196124 1 3
  3.16514 .9244325000000001 1 5
  3.19014 .9078579999999999 0 2
  3.22374 .9252921000000001 1 3
  3.28495          .9173229 1 3
  3.42511          .8885331 0 3
  3.47633 .8793797000000001 0 6
1.6195818 .9830000305175781 1 4
 1.629391               .98 1 4
  3.60902 .8759551000000001 1 4
   3.6502          .8555083 0 3
  3.66973 .9015787000000001 1 4
   1.6749 .9869999694824219 1 2
1.6876818 .9840000152587891 1 3
  3.71767          .8732349 0 3
  3.73397          .8507866 1 3
  3.77264           .864858 1 2
1.7370318 .9840000152587891 0 3
  3.82153 .8926086999999999 1 3
 1.749059              .985 1 2
  3.85073            .87676 0 3
  3.86381 .9192615000000001 1 1
  3.88459          .8217481 1 4
1.7716454 .9780000305175781 0 4
  3.91496 .8649952000000001 0 2
   3.9438          .8733145 0 2
  3.96859          .8162379 1 5
 1.811491 .9830000305175781 1 4
  3.98909           .874633 0 2
  4.00001          .8137358 1 4
  4.00957          .8648334 0 3
  4.01221 .8458074000000001 1 3
  4.01406          .8849974 0 2
  4.02415          .8647704 0 2
  4.04549 .8642723999999999 0 2
  4.05009          .8692121 1 2
 1.859109 .9790000152587891 1 4
   4.1098           .887962 0 2
 1.873932 .9830000305175781 1 2
  4.13317          .8747856 1 2
1.8896636 .9719999694824218 1 3
1.9208592 .9790000152587891 1 2
 1.922491               .98 0 2
  4.23509          .8130789 0 2
  4.24629 .8429949000000001 1 5
  4.25622          .8000713 0 4
   4.2659          .8567741 1 5
1.9608545 .9780000305175781 1 2
  4.32399 .8504305999999999 0 2
  4.39443 .7854139000000001 1 2
  4.40424          .8844606 1 2
2.0023546               .98 0 3
2.0034635  .975999984741211 0 4
  4.41624 .8245646999999999 0 2
  4.43171 .8841536000000001 1 4
2.0157363  .975999984741211 1 3
  4.45371          .8117735 1 3
   2.0397  .975999984741211 1 3
  4.51783          .7971114 0 2
2.0566409 .9769999694824219 0 3
 2.057718 .9719999694824218 1 4
 2.061068  .970999984741211 1 3
  4.54413 .8513078000000001 1 2
end
label values sex sex
label def sex 0 "female", modify
label def sex 1 "male", modify

[Joseph Coveney]

With the implied intercept (constant) that the main effects of sex term imposes, your nonlinear functional specification no longer fits the data (b0 is no longer different from zero).

So, if your nonlinear model isn't supposed to have a constant (most that are fitted with nl don't), then probably you should omit the main effects of sex term (the {s} * sex term) from the model specification that I suggested.

Your full dataset (N = 210) isn't included.

If you want to do that, I recommend not using dataex, which I believe is just to show a smattering of observations so that readers can get a feel for the thing. You should just attach the Stata dataset, perhaps as a comma-delimited or tab-delimited text file (help export delimited) for those Statalist members who are wary of viruses in binary-file attachments.

[Peta Hitchens]

Thanks again Joseph, greatly appreciated! The curves and fits are looking pretty good now. Problem is solved

Thanks also for the data example tip, I didn't notice it was only the first 100 using dataex, argh.
Have a lovely weekend,
/Peta