# Fit model with two dependent variables

Say I have an experiment that produces two values, for a given input. For example you get a complex response value.

Now I have a model that - of course - also produces two dependent values for one independent value, and some other parameters and I want to find the set of parameters that describe this model best (least-squares).

I basically need two equations for this model then (simple example)

$$Y_1(x) = A x + B$$ $$Y_2(x) = (A / 2) x + B$$

I want to find A and B, such that the overall sum of squares is minimal

$$\sum_i[(Y_{1,i} - Y_1(x_i))^2 + (Y_{2,i} - Y_2(x_i))^2] \Rightarrow \textrm{min}$$

and I also want to calculate the $$R^2$$ value for this fit.

How can I do this in Mathematica? I tried using the NonlinearModelFit function like this (data is in the form {X, Y1, Y2}):

data = {{1, 2, 1}, {2, 6, 7}, {3, 8, 9}, {4, 6, 5}, {5, 10, 12}, {6,
23, 18}};
NonlinearModelFit[data, {A*x + B, (A/2)*x + B}, {A, B}, x]


However I get the (not completely unexpected) error

Number of coordinates (2) is not equal to the number of variables (1).


The result I got from another program is A = 2.70824951933389 and B = 1.80751165636783 with R² = 0.52.

I just tried supplying two model functions instead of one. I couldn't find any examples for my case either, since you mostly find multiple-regression problems, with more than one independent variable instead, like in this question How can I use nonlinear model fit with multiple variables? .

• How is the correlation coefficient defined in that "another program"? (And what "another program" are you using, BTW?) Jul 10, 2015 at 17:00
• @belisarius I used Origin. From what I found by calculating it by hand they calculate the coefficient by dividing the sum of squared residuals of the fit by the sum of squared differences of each Y point (both Y1 and Y2) from the overall Y-mean (over all Y1 and Y2 values). I don't know if that is the correct approach though and I want to cross-check the method. Please see my other question on stats-overflow stats.stackexchange.com/questions/160326/…
– Jens
Jul 10, 2015 at 17:07
• For complex data specifically, I wrote a package to do the process explained in the answers automatically. Jul 10, 2015 at 20:55
• @OleksandrR. Nice work! I will give it a try
– Jens
Jul 11, 2015 at 6:07

There is no built-in way to accomplish what you want with NonlinearModelFit, or with the other fitting functions, unfortunately. Of course, you could write your own target function and NMinimize that, as belisarius has mentioned already in his answer. This is very general and not particularly difficult to do, but it does not give you access to the wealth of built-in fit description parameters that are accessible through NonlinearModelFit. Again, you could recalculate those manually, but it would be useful to be able to leverage the existing functionality instead.

In order to "trick" NonlinearModelFit into fitting your multivariate model, you will want to build an object function that includes both models and can "switch between them", based on a dummy selector variable that you will introduce in your data set and treat as an independent variable.

Let's go through a simple example. I generated some noisy data as a function of a common $x$ range which is the actual independent variable in your physical model. The two functions I used to generate the dependent variables are $y_1=x^2-5x+3$ and $y_2=0.1\ x^3-3$. I then added some random noise to the $y$ variables just to make this slightly more realistic. Here is the data set:

originaldata = {{-10, 168, -83}, {-9, 150.1, -86.9}, {-8, 129.6, -36.2}, {-7, 101.5, -34.3}, {-6, 74.8, -14.6}, {-5, 52.5, -32.5}, {-4, 48.6, -12.4}, {-3, 21.1, -23.7}, {-2, 22, -7.8}, {-1, 15.3, 5.9}, {0, 5, 8}, {1, -3.9, -7.9}, {2, -11.4, -7.2}, {3, -11.5, -11.3}, {4, 10.8, 22.4}, {5, 4.5, 8.5}, {6, 14.6, 29.6}, {7, 28.1, 31.3}, {8, 33, 31.2}, {9, 61.3, 68.9}, {10, 76, 96}};

ListPlot[originaldata[[All, {1, #}]] & /@ {2, 3}, PlotRange -> All]


You will have to reorganize this data set so that each point has the following form: {xValue, datasetIndentifier, yValue}, where we will use the numbers $1$ or $2$ as dataset indentifiers.

Transpose@originaldata;
Insert[%[[{1, #}]], ConstantArray[# - 1, Length[originaldata]], 2] & /@ {2, 3};
reorganizeddata = Flatten[Transpose /@ %, 1]


{{-10, 1, 153}, {-9, 1, 131}, {-8, 1, 107}, {-7, 1, 80}, {-6, 1, 59}, {-5, 1, 59}, {-4, 1, 31}, {-3, 1, 27}, {-2, 1, 21}, {-1, 1, 14}, {0, 1, 5}, {1, 1, 9}, {2, 1, -9}, {3, 1, -10}, {4, 1, 7}, {5, 1, 2}, {6, 1, 8}, {7, 1, 25}, {8, 1, 26}, {9, 1, 33}, {10, 1, 45}, {-10, 2, 139.36}, {-9, 2, 79.6883}, {-8, 2, 0.47189}, {-7, 2, -35.8288}, {-6, 2, -31.6333}, {-5, 2, -47.4649}, {-4, 2, -15.9864}, {-3, 2, 3.2233}, {-2, 2, -4.73885}, {-1, 2, -3.77884}, {0, 2, -20}, {1, 2, -0.221163}, {2, 2, -13.2612}, {3, 2, 13.7767}, {4, 2, 24.9864}, {5, 2, 52.4649}, {6, 2, 54.6333}, {7, 2, 12.8288}, {8, 2, -6.47189}, {9, 2, -80.6883}, {10, 2, -154.36}}

Now let's create an appropriate model function.

Clear[modelfunction]
modelfunction[indepvar_?NumericQ, datasetselector_?NumericQ, p11_?NumericQ, p12_?NumericQ, p13_?NumericQ, p21_?NumericQ, p22_?NumericQ] :=
Piecewise[{
{p11 indepvar^2 + p12 indepvar + p13,
datasetselector == 1},
{p21 indepvar^3 + p22,
datasetselector == 2}
}]


... and use it with NonlinearModelFit:

nlm = NonlinearModelFit[
reorganizeddata,
modelfunction[x, selector, a, b, c, d, e],
{a, b, c, d, e},
{x, selector}
]


We can then very conveniently extract any of the fit descriptors available through NonlinearModelFit:

nlm["ParameterTable"]
nlm["RSquared"]


Of course you can also plot the experimental points and the fit:

Plot[
nlm["BestFit"] /. selector -> # & /@ {1, 2}, {x, -10, 10},
Epilog -> {PointSize[0.015], Gray, Point[reorganizeddata[[All, {1, 3}]]]},
Evaluated -> True
]


• That is a very nice way (goes in line with the answer to a different question i found a few minutes ago, mentioned in my comment to belisarius answer). What I find curious, is that the R-Squared value I get from the nlm result is a lot larger than what Origin gives me (0.84 for my example data, as compared to 0.51 in Origin). Is the value, that comes out of the nlm correct, or is it distorted somehow by the workaround?
– Jens
Jul 10, 2015 at 18:06
• I found the difference in R²: groups.google.com/forum/#!searchin/… Seems like it's a matter of definition and Mathematica actually has a more general description of R² for nonlinear models. How are statisticians not going completely insane after 3 days?
– Jens
Jul 10, 2015 at 18:28
• @Jens That's an interesting point! I can't seem to open the link you provided though: it leads to a banned content warning page from Google. Do you have another way to access the content? Also, yeah, I have come to the conclusion that statisticians have a higher pain tolerance than I do, or otherwise they enjoy the pain they inflict upon others that try to use statistics... :-) Jul 10, 2015 at 18:37
• goo.gl/eEdcbJ Let's try that again, I couldn't either, even though I copy pasted the link. Ok that seems to work. Stackexchange added some part's to the URL for whatever reason that caused it to not work.
– Jens
Jul 10, 2015 at 18:38
• @Jens Your most recent link works for me. That's indeed an interesting point about Rsquared. Jul 10, 2015 at 18:57

If you know the form you want to minimize, Mathematica provides several ways to do that. For example:

data = {{1, 2, 1}, {2, 6, 7}, {3, 8, 9}, {4, 6, 5}, {5, 10, 12}, {6, 23, 18}};
y1[A_, B_, x_] := A*x + B
y2[A_, B_, x_] := (A/2)*x + B
s = Sum[(i[[2]] - y1[A, B, i[[1]]])^2 + (i[[3]] - y2[A, B, i[[1]]])^2, {i, data}];
NMinimize[s, {A, B}]
(* {211.085, {A -> 2.70825, B -> 1.80751}} *)


Which seems what Origin is doing

Edit

You may try to define something like a Pearson's coefficient:

m = Mean /@ Transpose@data;
mx = First@m;
my = m[[2 ;; 3]];
rr = {Sum[(i[[1]] - mx) (i[[2 ;; 3]] - my), {i, data}],
Sqrt[(Sum[(i[[1]] - mx)^2, {i, data}] Sum[(i[[2 ;; 3]] - my)^2, {i, data}])]};
Divide @@ (Norm /@ rr) // N
(* 0.855367 *)

• Thank you for the answer. I just found a different question about the same topic here: mathematica.stackexchange.com/questions/15905/… Getting the fitted values works well in both ways (yours and their answers). However the R² value I get from the NonlinearModelFit when using @b.gatessucks answer is different than the Origin answer. Can you show me a way to calculate R² using your approach?
– Jens
Jul 10, 2015 at 17:53
• @Jens See edit, please Jul 10, 2015 at 18:17
• That's a good idea! Due to the additional benefits of the NonlinearModelFit results I think I will use that one, as explained in the other answer, though. As you have pointed out correctly, there seem to be quite a lot of ways to do this.
– Jens
Jul 10, 2015 at 18:31
• @Jens I believe the other way has its shortcomings too.Just think about a case where due to multiplication factors the y2 series is orders of magnitude larger than the other. Jul 10, 2015 at 18:34
• The error function to minimize should be the same in both approaches, isn't it?
– Jens
Jul 10, 2015 at 18:37

If you simply have two sets of data and no shared parameters, one should just perform two fits.

If you have two separate sets of data with shared parameters and know that the errors are independent with a single common variance, then using MultiNonlinearModelFit by @SjoerdSmit is what you want.

But you appear have "paired" response variables (or "matched" response variables when there are more than 2) which (to me) screams that the error terms are likely correlated and maybe of unequal variance. And you have common parameters. For such cases the structure of the error terms become much more important and ignoring that structure will get you wrong estimates of precision of the non-error terms. (It is true that many times the estimates of the non-error terms won't differ by much even if you ignore the error structure. But one should really want an appropriate set of estimates of precision.)

Below I describe a reasonable starting model from which one should then check on the assumptions I'm going to make explicit.

We assume that the error structure is that of a bivariate normal distribution with that same distribution for all observations and that individual observations are independent of each other. (Again, one should check on those assumptions after fitting.)

Suppose @MarcoB 's models are slightly modified:

f1[x_, a_, b_, c_] := a x^2 - b x + c
f2[x_, c_, d_] := d x^3 - c

(* Generate some data *)
x = Range[-10, 10];
parms = {a -> 1, b -> 5, c -> 3, d -> 0.1, σ1 -> 7, σ2 -> 15, ρ -> 0.5};
(* Random error from bivariate normal *)
SeedRandom[12345];
error = RandomVariate[BinormalDistribution[{0, 0}, {σ1, σ2}, ρ] /. parms, Length[x]];
(* Construct response vectors *)
response = Transpose[{f1[x, a, b, c], f2[x, c, d]} /. parms] + error ;

(* Look at the data *)
ListPlot[{Transpose[{x, response[[All, 1]]}], Transpose[{x, response[[All, 2]]}]}]


(* Get initial estimates of parameters by fitting separately and ignoring the
correlation structure *)
(* The example functions are linear in the parameters but we'll use
NonlinearModelFit to keep things more general *)
nlm1 = NonlinearModelFit[Transpose[{x, response[[All, 1]]}], f1[z, a, b, c], {a, b, c}, z];
nlm2 = NonlinearModelFit[Transpose[{x, response[[All, 2]]}], f2[z, c, d], {c, d}, z];
{a0, b0, c10} = {a, b, c} /. nlm1["BestFitParameters"];
σ10 = Sqrt[nlm1["EstimatedVariance"]];
{d0, c20} = {d, c} /. nlm2["BestFitParameters"];
c0 = (c10 + c20)/2;  (* Average the two separate estimates of c *)
σ20 = Sqrt[nlm2["EstimatedVariance"]];
ρ0 = Correlation[Transpose[{nlm1["FitResiduals"], nlm2["FitResiduals"]}]][[1, 2]];

(* Define the log of the likelihood *)
logL = LogLikelihood[BinormalDistribution[{0, 0}, {σ1, σ2}, ρ],
response - Transpose[{f1[x, a, b, c], f2[x, c, d]}]];

(* Find maximum likelihood estimates *)
mle = FindMaximum[{logL, σ1 > 0 && σ2 > 0 && -1 < ρ < 1},
{{a, a0}, {b, b0}, {c, c0}, {d, d0}, {σ1, σ10}, {σ2, σ20}, {ρ, ρ0}}];

(* Estimate parameter covariance matrix and associated standard errors *)
cov = -Inverse[(D[logL, {{a, b, c, d, σ1, σ2, ρ}, 2}]) /. mle[[2]]];
se = Sqrt[Diagonal[cov]] ;(* Standard errors *)

(* Create a parameter table *)
t = ({a, b, c, d, σ1, σ2, ρ} /. mle[[2]])/se;
pvalue = 2 (1 - CDF[NormalDistribution[0, 1], Abs[t]]);
TableForm[Transpose[{{a, b, c, d, σ1, σ2, ρ} /. mle[[2]], se, t, pvalue}],
TableHeadings -> {{"a", "b", "c", "d", "σ1", "σ2", "ρ"},
{"Estimate", "Standard Error", "z-Statistic", "P-value"}}]


We can compare this against the parameter tables of the separate fits:

nlm1["ParameterTable"]


nlm2["ParameterTable"]


We see that the standard error for $$c$$ is about half of that of the separate fits.

Here is the data and the fit:

Show[ListPlot[{Transpose[{x, response[[All, 1]]}], Transpose[{x, response[[All, 2]]}]}],
Plot[{f1[z, a, b, c] /. mle[[2]], f2[z, c, d] /. mle[[2]]}, {z, Min[x], Max[x]}]]


With this approach one can also fit relatively simple mixed models.