# How can I combine NLS models in a way as is done with linear regression models in Path Analysis

I have a structural equation model that includes several paths that are actually non-linear (built in with the lavaan package for R). I have used polynomials in order to model these non-linear paths. Due to some variance issues I have to transform the variables which causes the polynomials to be non-significant.

In order to model these non-linear relationships I have taken several parts from my structural equation model and turned them into three non-linear least square (nls) models.

Now I wish to combine these nls models in order to calculate the total effect of one of the variables (V1) on another variable (v4). One of the models is actually a linear model but I created it with the nls package, because it is an important link. In the figure below, you can see the 3 different models and the conceptual path model. p1 <- nls (sqrt(V2) ~ (a*(1+V1)^-b),
start = list(a= 1, b=1), trace=TRUE,  control = list(maxiter = 500))
p2 <- nls (sqrt(V3) ~ (a*(V2)^-b),
start = list(a= 0.01219, b=0.62064), trace=TRUE,  control = list(maxiter = 500))
p3 <- nls (log(1+V4) ~ a+b*V2+c*V3,
start = list(a=0, b=1, c=1), trace=TRUE,data=data)


Note: that the response variables have been transformed in order to account for the non-equal variance.

I now wish to calculate the effect of changes in V1 on V4. In order to do that I would like to turn these nls models in a structural equation model or a path model. But I am not sure if this can be done and how this should be done. Can I just use the predicted values of each path as an input for the next path or am I missing something? Any help would be greatly appreciated.

Below follows the data, in case you would like to give it a go.

V1 <- c(0,0.0032,0.00414,0,0.00231,0.00408,0.00691,0,0,0.0026,0,0,0.00456,0.01071,0.01929,0,0,0.00351,0.0038,0.00693,0.00829,0.00714,0.01163,0.01389,0,0.01893,0,0.00586,0,0.0037,0.00598,0.02007,0.03342,0,0,0,0.00291,0.0042,0.0069,0,0.00322,0.0332,0,0,0,0,0,0,0.00554,0,0,0,0.0025,0.00521,0.00659,0.01577,0.02058)
V2 <- c(0.0521094298,0.0447902107,0.0266698315,0.0604551725,0.0184757506,0.0225835592,0.0207404335,0.0626592589,0.1457083799,0.0519345625,0.1267201081,0.1218615401,0.0119282397,0.0214224507,0.0192882631,0.2739833262,0.0666306861,0.0281081461,0.0972920383,0.0069295507,0.0538835535,0.0136833223,0.0087190369,0.0104191267,0.2157022429,0.0283973259,0.0595832837,0.0740827387,0.0730018778,0.0296174536,0.0475871169,0.0267537054,0.0200494553,0.1351056422,0.0894515929,0.0881266729,0.0267599698,0.0117233294,0.0620573244,0.0771381738,0.0177237201,0.0664084117,0.0777634875,0.3840885868,0.238855232,0.14546571,0.1240576391,0.093863663,0.12182966,0.3828426666,0.3512659359,0.0696672472,0.105055154,0.0399194668,0.1021624383,0.0315415955,0.0051445357)
V3 <- c(0.044,0.038,0.033,0.022,0.231,0.192,0.062,0.224,0.034,0.039,0.009,0.004,0.115,0.161,0.386,0.027,0.087,0.246,0.011,0.159,0.033,0.579,0.084,0.049,0.022,0.062,0.055,0.042,0.049,0.259,0.208,0.054,0.12,0.016,0.007,0.002,0.131,0.263,0.159,0.12,0.006,0.066,0.026,0,0,0.028,0.052,0.012,0,0.072,0.08,0.082,0.1,0.05,0.011,0.226,0.412)
V4 <- c(10,9,6.6666666667,6,5.3333333333,5,4,3.5,3,3,2.6666666667,2.6666666667,2.6666666667,2.6666666667,2.6666666667,2.5,2.5,2.3333333333,2.3333333333,2.3333333333,2.3333333333,2,2,2,1.6666666667,1.6666666667,1.5,1.5,1.3333333333,1.3333333333,1.3333333333,1.3333333333,1.3333333333,1,1,1,1,1,1,0.6666666667,0.6666666667,0.6666666667,0.5,0.3333333333,0.3333333333,0.3333333333,0.3333333333,0.3333333333,0.3333333333,0,0,0,0,0,0,0,0)


Your estimated model looks like: \begin{align} \sqrt{V_2} &= a_2 \left( 1+V_1 \right)^{b_2} + \epsilon_2 \\ \sqrt{V_3} &= a_3 \left( 1+V_2 \right)^{b_3} + \epsilon_3 \\ \ln{\left(1+V_4\right)} &= a_4 + b_4 V_2 + c_4 V_3 + \epsilon_4 \\ \end{align}

What you are interested in calculating is $\frac{d V_4}{d V_1}$. To do this, we are just going to use the Chain Rule: \begin{align} \frac{d V_4}{d V_1} &= \frac{\partial V_4}{\partial V_2}\frac{d V_2}{d V_1}+\frac{\partial V_4}{\partial V_3}\frac{d V_3}{d V_2}\frac{d V_2}{d V_1}\\ \end{align} To calculate those derivatives, we have to re-write the equations you estimated so that they have the correct variables on the left-hand-side: \begin{align} V_2 &= \left( a_2 \left( 1+V_1 \right)^{b_2} + \epsilon_2 \right) ^2\\ V_3 &= \left( a_3 \left( 1+V_2 \right)^{b_3} + \epsilon_3 \right) ^2\\ V_4 &= \exp{ \left(a_4 + b_4 V_2 + c_4 V_3 + \epsilon_4 \right)} -1\\ \end{align} Now, we take the derivatives: \begin{align} \frac{d V_2}{d V_1} &= 2\left( a_2 \left( 1+V_1 \right)^{b_2} + \epsilon_2 \right) \left( a_2 b_2 \left( 1+V_1 \right)^{b_2-1} \right)\\ &\\ \frac{d V_3}{d V_2} &= 2 \left( a_3 \left( 1+V_2 \right)^{b_3} + \epsilon_3 \right) \left( a_3 b_3 \left( 1+V_2 \right)^{b_3-1} \right)\\ &\\ \frac{\partial V_4}{\partial V_2} &= \exp{ \left(a_4 + b_4 V_2 + c_4 V_3 + \epsilon_4 \right)}b_4\\ &\\ \frac{\partial V_4}{\partial V_3} &= \exp{ \left(a_4 + b_4 V_2 + c_4 V_3 + \epsilon_4 \right)}c_4\\ \end{align} Now, we can use the original equations to make some substitutions: \begin{align} \frac{d V_2}{d V_1} &= 2 \sqrt{V_2} \left( a_2 b_2 \left( 1+V_1 \right)^{b_2-1} \right)\\ &\\ \frac{d V_3}{d V_2} &= 2 \sqrt{V_3} \left( a_3 b_3 \left( 1+V_2 \right)^{b_3-1} \right)\\ &\\ \frac{\partial V_4}{\partial V_2} &= b_4\ln{\left(1+V_4\right)}\\ &\\ \frac{\partial V_4}{\partial V_3} &= c_4\ln{\left(1+V_4\right)} \end{align} Obviously, those have $V_1,V_2,V_3,V_4$ in them. This is as you should expect---in non-linear models, the derivative (the effect size) depends on the levels of the variables. That's what non-linear means.

Since the derivatives depend on data, the value of the derivative will be different for each data point. Again, as expected. Now, since you usually want to report one number for your effects, the number you usually report is the average of the effect sizes. So, what you do is calculate the following for each data point: \begin{align} \frac{d V_4}{d V_1} &= \frac{\partial V_4}{\partial V_2}\frac{d V_2}{d V_1}+\frac{\partial V_4}{\partial V_3}\frac{d V_3}{d V_2}\frac{d V_2}{d V_1}\\ &\\ \frac{d V_4}{d V_1} &= b_4\ln{\left(1+V_4\right)} \left(2 \sqrt{V_2} \left( a_2 b_2 \left( 1+V_1 \right)^{b_2-1} \right) \right)\\ &+ c_4\ln{\left(1+V_4\right)} \left( 2 \sqrt{V_3} \left( a_3 b_3 \left( 1+V_2 \right)^{b_3-1} \right) \right) \left(2 \sqrt{V_2} \left( a_2 b_2 \left( 1+V_1 \right)^{b_2-1} \right) \right)\\ & \\ \frac{d V_4}{d V_1} &= 2 b_4 a_2 b_2\ln{\left(1+V_4\right)} \sqrt{V_2} \left( 1+V_1 \right)^{b_2-1}\\ &+ 4 c_4 a_3 b_3 a_2 b_2 \ln{\left(1+V_4\right)} \sqrt{V_3} \left( 1+V_2 \right)^{b_3-1} \sqrt{V_2} \left( 1+V_1 \right)^{b_2-1} \end{align} The last equation is what you use to calculate, for each observation, the total effect of $V_1$ on $V_4$. The first big term (before the $+$ sign) is the effect of $V_1$ on $V_4$ which operates through $V_2$ directly, and the second big term (after the $+$ sign) is the effect of $V_1$ on $V_4$ which operates through $V_2$ and then $V_3$.

To present a single number for the overall effect, you just take the sample mean over the observations in your dataset. Similarly, for the effects running through $V_2$ and $V_3$, you take sample means over their respective terms.

Now, you probably will want to construct a confidence interval for these effects as well. The easiest and probably best way to do this is via bootstrapping. The R code below does all of this.

Notice that there is a fair amount of algebra above and a fair amount of translating algebra into code necessary to run the R code. It would not be shocking if there were errors.

On the bright side, the median effect size my code calculates for the total effect of V1 on V4 is -38. A simple linear regression of V4 on V1 gives -32. So, that's kind of reassuring.

# This script written in response to a Cross Validated question
# http://stats.stackexchange.com/questions/264259/how-can-i-combine-nls-models-in-a-way-as-is-done-with-linear-regression-models-i

# The script estimates a three-equation structural equation model using OLS
# and NLS.  Then, it calculates the total effect of variable v1 on variable
# v4.  This total effect operates (non-linearly) through two intermediary
# variables, v2 and v3.  Then, the script calculates a confidence
# interval via bootstrapping.

# William B Vogt

library(Hmisc)
library(boot)
set.seed(12344321)

# data from the question
V1 <- c(0,0.0032,0.00414,0,0.00231,0.00408,0.00691,0,0,0.0026,0,0,0.00456,0.01071,0.01929,0,0,0.00351,0.0038,0.00693,0.00829,0.00714,0.01163,0.01389,0,0.01893,0,0.00586,0,0.0037,0.00598,0.02007,0.03342,0,0,0,0.00291,0.0042,0.0069,0,0.00322,0.0332,0,0,0,0,0,0,0.00554,0,0,0,0.0025,0.00521,0.00659,0.01577,0.02058)
V2 <- c(0.0521094298,0.0447902107,0.0266698315,0.0604551725,0.0184757506,0.0225835592,0.0207404335,0.0626592589,0.1457083799,0.0519345625,0.1267201081,0.1218615401,0.0119282397,0.0214224507,0.0192882631,0.2739833262,0.0666306861,0.0281081461,0.0972920383,0.0069295507,0.0538835535,0.0136833223,0.0087190369,0.0104191267,0.2157022429,0.0283973259,0.0595832837,0.0740827387,0.0730018778,0.0296174536,0.0475871169,0.0267537054,0.0200494553,0.1351056422,0.0894515929,0.0881266729,0.0267599698,0.0117233294,0.0620573244,0.0771381738,0.0177237201,0.0664084117,0.0777634875,0.3840885868,0.238855232,0.14546571,0.1240576391,0.093863663,0.12182966,0.3828426666,0.3512659359,0.0696672472,0.105055154,0.0399194668,0.1021624383,0.0315415955,0.0051445357)
V3 <- c(0.044,0.038,0.033,0.022,0.231,0.192,0.062,0.224,0.034,0.039,0.009,0.004,0.115,0.161,0.386,0.027,0.087,0.246,0.011,0.159,0.033,0.579,0.084,0.049,0.022,0.062,0.055,0.042,0.049,0.259,0.208,0.054,0.12,0.016,0.007,0.002,0.131,0.263,0.159,0.12,0.006,0.066,0.026,0,0,0.028,0.052,0.012,0,0.072,0.08,0.082,0.1,0.05,0.011,0.226,0.412)
V4 <- c(10,9,6.6666666667,6,5.3333333333,5,4,3.5,3,3,2.6666666667,2.6666666667,2.6666666667,2.6666666667,2.6666666667,2.5,2.5,2.3333333333,2.3333333333,2.3333333333,2.3333333333,2,2,2,1.6666666667,1.6666666667,1.5,1.5,1.3333333333,1.3333333333,1.3333333333,1.3333333333,1.3333333333,1,1,1,1,1,1,0.6666666667,0.6666666667,0.6666666667,0.5,0.3333333333,0.3333333333,0.3333333333,0.3333333333,0.3333333333,0.3333333333,0,0,0,0,0,0,0,0)

# estimation technique from the question (slightly edited)
p2 <- nls (sqrt(V2) ~ (a*(1+V1)^-b),
start = list(a= 1, b=1), trace=TRUE,  control = list(maxiter = 500))
p3 <- nls (sqrt(V3) ~ (a*(V2)^-b),
start = list(a= 0.01219, b=0.62064), trace=TRUE,  control = list(maxiter = 500))
p4 <- lm (log(1+V4) ~ V2 + V3)

# look at coefficients
p2
p3
p4

# Now calculate total effect, dV4dV1 and the two partial effects, direct2 and thru3
direct2 <- 2*coef(p4)["V2"]*coef(p2)["a"]*coef(p2)["b"]*log(1+V4)*sqrt(V2)*(1+V1)^(coef(p2)["b"]-1)
thru3 <- 2*coef(p4)["V3"]*coef(p3)["a"]*coef(p3)["b"]*log(1+V4)*sqrt(V3)*(1+V2)^(coef(p3)["b"]-1)*direct2
dv4dv1 <- direct2 + thru3

# OK, let's look at our effects.  Describe will provide us with mean, median, and several
# percentiles of the effect
describe(cbind(direct2,thru3,dv4dv1))

# OK, now let's bootstrap in order to get confidence intervals
# This requires writing a routine which calcuates and returns the three estimators above
# To do this, it's convenient to put the data into a data frame

my.data <- data.frame(V1,V2,V3,V4)

# The function to do one bootstrap replication
effect.calc <- function(data,eq2,eq3,indices){
# Ensure args supplied
if(missing(data)){stop("No data supplied")}
if(missing(eq2)){stop("No eq2 supplied")}
if(missing(eq3)){stop("No eq3 supplied")}

# Create bootstrap sample
boot.data <- data[indices,]

# Run estimations, using converged values as starting values
p2.bs <- nls (sqrt(V2) ~ (a*(1+V1)^-b),
start = coef(eq2), trace=FALSE,  control = list(maxiter = 500),data=boot.data)
p3.bs <- nls (sqrt(V3) ~ (a*(V2)^-b),
start = coef(eq3), trace=FALSE,  control = list(maxiter = 500),data=boot.data)
p4.bs <- lm (log(1+V4) ~ V2+V3,data=boot.data)

# Calculate effects
direct2 <- 2*coef(p4.bs)["V2"]*coef(p2.bs)["a"]*coef(p2.bs)["b"]*log(1+boot.data$V4)*sqrt(boot.data$V2)*(1+V1)^(coef(p2.bs)["b"]-1)
thru3 <- 2*coef(p4.bs)["V3"]*coef(p3.bs)["a"]*coef(p3.bs)["b"]*log(1+boot.data$V4)*sqrt(boot.data$V3)*(1+V2)^(coef(p3.bs)["b"]-1)*direct2
dv4dv1 <- direct2 + thru3

# Return mean effects
return(c(effect=mean(dv4dv1),direct=mean(direct2),thru=mean(thru3)))

}

# Do 10,000 bootstrap replications and store in dataframe boots
boots <- boot(data=my.data,statistic=effect.calc,R=10000,eq2=p2,eq3=p3)

boots

boot.ci(boots, type = "all", index = 1)  # overall effect of v1 on v4
boot.ci(boots, type = "all", index = 2)  # direct through v2
boot.ci(boots, type = "all", index = 3)  # indirect through v3

• @ Bill Thanks for your elaborate answer. It has been quite helpful, but I am still a bit figuring out all the details about what you answered. Getting a single measure for the effects is one thing, now I am trying to figure out how to use this properly to predict V4 based on changes in V1. I also noticed that I am a bit late with replying and therefore I am not sure what happened to the bounty. Did you receive it? – Robbie Mar 14 '17 at 6:05
• If you mean predict changes in $V_4$ based on changes in $V_1$, then the formulas and code I gave above do that. The mean effect I calculate gives the average amount (over all the data points) that $V_4$ goes up for a one unit increase in $V_1$. If there is a particular data point for which you are interested in calculating the effect on $V_4$ of $V_1$ going up by one, then just look at dv4dv1 for that point. If you need to predict the level of $V_4$ given the level of $V_1$, let me know. I didn't receive the bounty, but it's not that big a deal. – Bill Mar 14 '17 at 12:38