# Multivariate multiple nonlinear regression in R

I want to do multivariate (with more than 1 response variables) multiple (with more than 1 predictor variables) nonlinear regression in R.

The data I am concerned with are 3D-coordinates, thus they interact with each other, i.e. the x,y,z-coordinates are not independent. So I cannot just call the nls separately for each response variable (which I tried at first).

A subset of the data-frame with 3D-coordinates where x,y,z are the predictive variables and a,b,c the response variables:

              x           y         z           a            b         c
1  -2.26470e-03 -0.05081670 0.0811701 -0.00671079 -0.045721600 0.0705679
2  -9.13106e-05 -0.00670734 0.0724838 -0.00676299 -0.001638430 0.0588486
3   3.81399e-04  0.03556000 0.0782059 -0.00783726  0.038503800 0.0641364
4   1.42293e-03  0.06133920 0.0708688 -0.00820760  0.062697100 0.0572740
5  -5.06043e-02  0.04759040 0.0418189 -0.05949350  0.040427800 0.0266159
6   5.92963e-02  0.04183450 0.0431029  0.05124780  0.038396500 0.0327903
7  -4.44213e-02 -0.00909717 0.0459059 -0.05021130 -0.005634520 0.0329833
8  -3.75400e-02 -0.00625770 0.0567296 -0.04255200 -0.000666089 0.0436465
9  -2.37768e-02 -0.00707318 0.0581552 -0.03048950 -0.001260670 0.0457355
10 -1.56645e-02 -0.01326670 0.0540247 -0.02101350 -0.009021990 0.0413755


My question: Is it possible to call the nls function with more than 1 (in my case 3) response variables? In other words is it possible to substitute y in nls(y ~ f(x,y,z, parameters), data) with something like c(a,b,c) or cbind(a,b,c), such that nls(cbind(a,b,c) ~ f(x,y,z, parameters), data) ?

In the post How to write R formula for multivariate response? it is shown that one can combine several response variables with cbind in the case of linear modeling with the lm function. This doesn't seem to work for nonlinear modeling with nls .., because the nls call in the code sample at the bottom of my question throws the following error:

Error in parse(text = x) : <text>:2:0: unexpected end of input 1: ~ ^

which I could not find a solution for online concerning my case of a multivariate regression..

My web-searches to my main question only gave me results concerning multivariate linear regression, which for example included solutions with the manova function..

Therefore, my question asked in a more general way: How do you in general solve such a non-linear multivariate multiple regression problem in R which takes into account interactions/dependencies between variables?

Here is my code where

• function f computes the rotations of coordinates about three axes in the order x-axis, y-axis, and then z-axis (unfortunately I cannot include the pic of the equation I wrote in LaTeX here since I haven't got 10 reputation points yet);
• rot_data_all is structured as the data-subset above, just with more rows;
• alpha1, alpha2 and so on are the parameters which nonlinear regression should approximate:

The code:

f <-function(x, y, z, alpha1, alpha2, alpha3, gamma, theta, phi, s) {
a <- alpha1 + s*(cos(theta)*cos(phi)*x - cos(theta)*sin(phi)*y + sin(theta)*z)
b <- alpha2 + s*((sin(gamma)*sin(theta)*cos(phi) + cos(gamma)*sin(phi))*x
+ (-sin(gamma)*sin(theta)*sin(phi) + cos(gamma)*cos(phi))*y
- sin(gamma)*cos(phi)*z)
c <- alpha3 + s*((cos(gamma)*sin(theta)*cos(phi) + sin(gamma)*sin(phi))*x
+ (cos(gamma)*sin(theta)*sin(phi) + sin(gamma)*cos(phi))*y
+ cos(gamma)*cos(phi)*z)
return(c(a,b,c))
}

rot.nls <- nls(cbind(a, b, c) ~ f(x, y, z, alpha1, alpha2, alpha3, gamma, theta, phi, s),
data = rot_data_all,
start = c(alpha1 = 0, alpha2 = 0, alpha3 = 0, gamma = 0.1, theta = 0.1, phi = 0.1, s = 0.1), trace = TRUE)


I hope to find a solution which is general enough to also solve other transformations which cannot be easily linearized like the set of equations for projective transformation, i.e. something like the following function:

f.proj <-function(x, y, z, betas) {
a <- (betas[1,1]*x + betas[1,2]*y + betas[1,3]*z + betas[1,4]) / (betas[4,1]*x + betas[4,2]*y + betas[4,3]*z + betas[4,4])
b <- (betas[2,1]*x + betas[2,2]*y + betas[2,3]*z + betas[2,4]) / (betas[4,1]*x + betas[4,2]*y + betas[4,3]*z + betas[4,4])
c <- (betas[3,1]*x + betas[3,2]*y + betas[3,3]*z + betas[3,4]) / (betas[4,1]*x + betas[4,2]*y + betas[4,3]*z + betas[4,4])
return(c(a,b,c))
}


I am happy to provide more information if needed! Thank you so much!

• You might look into the nice book istics.net/pdfs/multivariate.pdf It has an appendix containing R code that is simple to implement. Mar 18 '15 at 15:52

A bit late, but for posterity, here is a solution using tidy (or skinny) data. Each coordinate variable can be thought of as a measurement of the same kind made in different dimensions.

con <- textConnection(
"              x           y         z           a            b         c
1  -2.26470e-03 -0.05081670 0.0811701 -0.00671079 -0.045721600 0.0705679
2  -9.13106e-05 -0.00670734 0.0724838 -0.00676299 -0.001638430 0.0588486
3   3.81399e-04  0.03556000 0.0782059 -0.00783726  0.038503800 0.0641364
4   1.42293e-03  0.06133920 0.0708688 -0.00820760  0.062697100 0.0572740
5  -5.06043e-02  0.04759040 0.0418189 -0.05949350  0.040427800 0.0266159
6   5.92963e-02  0.04183450 0.0431029  0.05124780  0.038396500 0.0327903
7  -4.44213e-02 -0.00909717 0.0459059 -0.05021130 -0.005634520 0.0329833
8  -3.75400e-02 -0.00625770 0.0567296 -0.04255200 -0.000666089 0.0436465
9  -2.37768e-02 -0.00707318 0.0581552 -0.03048950 -0.001260670 0.0457355
10 -1.56645e-02 -0.01326670 0.0540247 -0.02101350 -0.009021990 0.0413755
")
# reshape data
library(tidyr)
rot_data_skinny <- gather(rot_data_all, key = coord, value = value, x, y, z)
rot_data_skinny <- gather(rot_data_skinny, key = coorda, value = valuea, a, b, c)
#   coord        value coorda      valuea
# 1     x -2.26470e-03      a -0.00671079
# 2     x -9.13106e-05      a -0.00676299
# 3     x  3.81399e-04      a -0.00783726
# 4     x  1.42293e-03      a -0.00820760
# 5     x -5.06043e-02      a -0.05949350
# 6     x  5.92963e-02      a  0.05124780


The function just needs a small tweak to accept the single input column.

f <- function(value, coord, alpha1, alpha2, alpha3, gamma, theta, phi, s) {
x <- value[coord == "x"]
y <- value[coord == "y"]
z <- value[coord == "z"]
a <- alpha1 + s*(cos(theta)*cos(phi)*x - cos(theta)*sin(phi)*y + sin(theta)*z)
b <- alpha2 + s*((sin(gamma)*sin(theta)*cos(phi) + cos(gamma)*sin(phi))*x +
(-sin(gamma)*sin(theta)*sin(phi) + cos(gamma)*cos(phi))*y -
sin(gamma)*cos(phi)*z)
c <- alpha3 + s*((cos(gamma)*sin(theta)*cos(phi) + sin(gamma)*sin(phi))*x +
(cos(gamma)*sin(theta)*sin(phi) + sin(gamma)*cos(phi))*y +
cos(gamma)*cos(phi)*z)
return(c(a,b,c))
}


The response is now just the position in space grouped by dimension. Note that the parameter significance estimates are probably being overestimated as the three measurements for position are being treated as three independent measurements. But in terms of being able to estimate a transform, this seems like it would have been useful...

rot.nls <- nls(valuea ~ f(value, coord, alpha1, alpha2, alpha3,
gamma, theta, phi, s),
data = rot_data_skinny,
start = c(alpha1 = 0, alpha2 = 0, alpha3 = 0, gamma = 0.1,
theta = 0.1, phi = 0.1, s = 0.1))
summary(rot.nls)
#
# Formula: valuea ~ f(value, coord, alpha1, alpha2, alpha3, gamma, theta,
#     phi, s)
#
# Parameters:
#          Estimate Std. Error t value Pr(>|t|)
# alpha1 -9.932e-03  8.134e-04 -12.211  < 2e-16 ***
# alpha2 -6.734e-05  7.509e-04  -0.090 0.928755
# alpha3 -8.769e-03  6.429e-04 -13.640  < 2e-16 ***
# gamma  -5.242e-02  1.114e-02  -4.705    1e-05 ***
# theta   4.564e-02  1.227e-02   3.719 0.000362 ***
# phi     1.597e-02  8.750e-03   1.826 0.071503 .
# s       9.493e-01  7.945e-03 119.478  < 2e-16 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.002041 on 83 degrees of freedom
#
# Number of iterations to convergence: 7
# Achieved convergence tolerance: 1.208e-06