# Independence of residuals in a computer-based experiment/simulation?

I conducted a computer-based assessment of different methods of fitting a particular type of model used in the palaeo sciences. I had a large-ish training set and so I randomly (stratified random sampling) set aside a test set. I fitted $m$ different methods to the training set samples and using the $m$ resulting models I predicted the response for the test set samples and computed a RMSEP over the samples in the test set. This is a single run.

I then repeated this process a large number of times, each time I chose a different training set by randomly sampling a new test set.

Having done this I want to investigate if any of the $m$ methods has better or worse RMSEP performance. I also would like to do multiple comparisons of the pair-wise methods.

My approach has been to fit a linear mixed effects (LME) model, with a single random effect for Run. I used lmer() from the lme4 package to fit my model and functions from the multcomp package to perform the multiple comparisons. My model was essentially

lmer(RMSEP ~ method + (1 | Run), data = FOO)


where method is a factor indicating which method was used to generate the model predictions for the test set and Run is an indicator for each particular Run of my "experiment".

My question is in regard to the residuals of the LME. Given the single random effect for Run I am assuming that the RMSEP values for that run are correlated to some degree but are uncorrelated between runs, on the basis of the induced correlation the random effect affords.

Is this assumption of independence between runs valid? If not is there a way to account for this in the LME model or should I be looking to employ another type of statical analysis to answer my question?

-
Are the residuals conditional on the predicted random effects or unconditional and in the simulations are the predicted random effects constant or varying. Remember try to get a sense of this for the default simulation methods in LME4 and not being able to (but the project was cancelled before I sorted it out). – phaneron Oct 18 '12 at 18:53
Not sure I follow fully, but the various runs of draw training set -> fit models -> compute RMSEP are all done prior to the LME. The random effect is for run as each run will have a different intercept (RMSEP) as different combinations of test set samples are chosen, but this is constant within run. As for the conditional/unconditional bit, I'm not sure/clear what you mean. Thanks for you comment. – Gavin Simpson Oct 18 '12 at 19:02

## 2 Answers

You are essentially doing some form of cross-validation here for each of your m methods and would then like to see which method performed better. The results between runs will definitely be dependent, since they are based on the same data and you have overlap between your train/test sets. The question is whether this should matter when you come to compare the methods.

Let's say you would perform only one run, and would find that one method is better than the others. You would then ask yourself - is this simply due to the specific choice of test set? This is why you repeat your test for many different train/test sets. So, in order to determine that a method is better than other methods, you run many times and in each run compare it to the other methods (you have different options of looking at the error/rank/etc). Now, if you find that a method does better on most runs, the result is what it is. I am not sure it is helpful to give a p-value to this. Or, if you do want to give a p-value, ask yourself what is the background model here?

-
Thanks for your thoughts. I think your last lines sums up pretty much where I am now. In anticipation of this somewhat I have a follow-up where I ask about appropriate ways of analysing this type of data. I also like your point about "it is what it is"; that had been swirling at the edges of my thought process recently too. – Gavin Simpson Oct 19 '12 at 18:47
One issue I have with the "result is what it is" part is that the RMSEPs are quite variable from run to run. So on average one or two methods are better, but are they really better given the variability in the RMSEPs? Hence my trying an LME with random effect for Run. To modify that approach I would need to know who correlated each data set is. It would seem that any statistical test I do would need to be so modified. Hence I still struggle with how to interpret the means from the 50 Runs for each method & whether I can draw any conclusions...? – Gavin Simpson Oct 22 '12 at 11:01
The way I see it, evaluating your methods over all possible train/test set partitions of you data would have been the most comprehensive evaluation. Since this is impossible, you are estimating this with random runs. Let's say you could evaluate all train/test partitions - you would still be left with the question of how to decide which method is better. So this is more a matter of how you define what "good" is. Does it mean high average score? Or does it mean that in many runs one method gets a higher score than the others (personally I think this would be a better version)? – Bitwise Oct 22 '12 at 13:23

May not really understand what you have done but

for Run I am assuming that the RMSEP values for that run are correlated to some degree

Yes, that reflects how challenging the test set was in that run

but are uncorrelated between runs

No, given the way you have sampled the test sets some will be more overlapped than others (most definitely not independent replications)

You would somehow have to model the dependency based on the overlap or design the assessment so the runs are independent. I would read the stats literature on cross-validation ;-)

-
+1 Thanks for the Answer. Hmm, I see what you mean. The more similar the test sets are the more similar their RMSEP values will be. OK, put that way it is the same as if the data were spatially or temporarily correlated. The way I generate the training sets / test sets should mean that on average they are all as dissimilar from one another. I'm not sure what CV would get me here - and in a sense I am doing that anyway just via a resampling approach. Will probably ask another Q then on how to solve the real problem. – Gavin Simpson Oct 19 '12 at 12:47
I'll leave this open until the end of the bounty period to see if anyone else bites, but I do appreciate your thoughts here and will accept and award bounty if no other Answers are forthcoming. – Gavin Simpson Oct 19 '12 at 12:48