# multivariate mean?

[edit : because my question was ambiguous, I decided to rewrite it entirely, with some simplification but a lot more details on the experimental design]

Four independent 10m*10m plots each received sewage sludge from the same water treatment facility (this is precisely where pseudoplication occurs).

After a year of waiting, one "representative" 1kg soil sample was sampled from each plot.

Ten 1g-sub-samples were sampled from each soil sample, then suspended in 10ml water and agitated for 2 days in order to "equilibrate".

Only then, 1 (arbitrary) unit of labeled molecule was added into each soil suspension.

We want to follow the decrease of labeled molecule quantity as time goes on.

After 1 min of agitation, remaining labeled molecule in water of one soil suspension was measured destructively, giving a value "v" between 0 and 1. The value of w = 1-v was recorded.

Same after 2, 3, 4, 5, 6, 7, 8, 9 and 10 minutes of agitation.

So far, we have recorded 10 w values for each soil sample.

plot of w with time of agitation :

Because v decreases with time, w is increasing with time and can be efficiently modeled by the following model equation :

      w ~ a * time^b


This model was fit onto each of the 4 sets of 10 w values, giving four sets of (a,b) parameters.

My problem is : How to calculate / estimate a set of global (a,b) parameters ?

I think calculating average of each parameter (=> a_mean and b_mean) is not right. Fitting a model on all pooled 40 w values is no more right.

NB : As a side note, if I had one, and only one, parameter of interest (for example, a maximum value, or a mean) by soil sample, one good option would be to average those four values into one global average.

Thanks !

• Do you have any idea (optimally: a distribution) how the w-vectors from different plots are different from each other? Can you provide the data or a plot? – A. Donda Jun 11 '15 at 15:32
• @A.Donda plot added. w-vectors are not really different in practice. – Rodolphe Jun 11 '15 at 16:53

Your model assumption is $$w = a ~ t^b$$ with some small random variations.

Take the logarithm of that equation: $$\log w = \log a + b ~ \log t.$$ You now have a linear relationship between $\log w$ and $\log t$ with two parameters, $\log a$ and $b$. Determining the parameters is a standard linear regression problem, and if for simplicity we assume normally distributed errors with constant variance, the parameters can be estimated by ordinary least squares.

I don't know which software you use, but in Matlab you would do the following: Define the time points:

t = (1 : 10)';


Let's assume that your $w$-values are contained in a matrix W with 10 rows for the different time points and 4 columns for the different samples. Set up a design matrix where the first column consists of repetitions of $\log t$ and the second is constant 1

X = [repmat(log(t), size(W, 2), 1), ones(numel(W), 1)];


and construct a vector of logarithmized data

y = log(W(:));


Now you can perform the linear regression using Matlab's backslash operator:

parest = X \ y;


The estimates of the parameters are

b = parest(1);
a = exp(parest(2));

• From what I understand, you seem to explain how to estimate as many (a,b) as columns in the W matrix. However my question is about how to get a global (a,b) while taking into account that soil samples are pseudoreplicates. – Rodolphe Jun 11 '15 at 20:36
• @Rodolphe, no, read again. You get exactly one estimate of (a,b). Try it out for yourself. – A. Donda Jun 11 '15 at 22:28
• I just achieved to try it by myself with FreeMat (I don't have matlab actually, and I am not familiar at all with the syntax so it took me some time... I use R). Now I understand that your way of doing this estimates one set of parameter by pooling all data together and not taking into account the relationship between datas from the same soil. I am sorry but this is exactly what I said in my question to be wrong. Thanks for showing me a bit of matlab syntax however. – Rodolphe Jun 12 '15 at 9:00
• @Rodolphe, well it is "right" if you are willing to make the assumption that you have normally distributed errors in the logarithmized relation. If you think it is not right, you should explain why so that an answerer can deal with the underlying problem; just arbitrarily excluding a legitimate solution doesn't make sense. Btw., the standard free alternative to Matlab is Octave, and all my code should run with it unchanged. – A. Donda Jun 22 '15 at 10:04
• In your own answer you assume that the true parameters are different and normally distributed; also a legitimate approach. But in that case you could just as well estimate the parameters on each data set and then average, which you also declared "not right". You're not playing fair here! :) – A. Donda Jun 22 '15 at 10:06

Thanks to some simple R simulations, I came up with a satisfying answer to my own question.

The question asks what is the right way of estimating the true / global set of parameters of a given model, given the model has been fitted to N different datasets, yielding N vectors of parameters.

The answer is : calculate vector average from these N vectors.

How I came up with this answer that contradicts what I though at first is quite well demonstrated by a simulation (= backing up my answer with personal experience) :

Let's imagine we observe two variables X (independent variable), and Y (dependent variable). Assuming Y is linearly dependent on X, we can construct the following linear model :

 Y = a * X + b


Assume we know by construction that a_TRUE = 2 and b_TRUE = 10.

In the case of a and b being orthogonal to each other, one can draw a sample of N "a" values, and another sample of N "b" values.

 a_sample <- rnorm(10,a_TRUE,1)
b_sample <- rnorm(10,b_TRUE,2)


here is one possible plot of the true equation (black solid) and all 10 individual equations (dashed red) :

By construction, a_TRUE can be estimated by averaging all individual vectors.

Here (sorry, I forgot to set seed so it is not reproducible), b_estimate = 9.7 +- 0.7 (95% CI), and a = 1.5 +- 0.8

By analogy, in the nonlinear case, it seems OK to average over all 4 (a,b) vectors.

However, with only four values to average by parameter, it is likely that 95% confident intervals will be huge.

I don't know, thought, if this is still OK when the parameters are correlated.