# Comparing the within-subject variance between two groups of subjects

If each subject has repeated a test multiple times, I am able to compute the mean and variance of their performance. (I am assuming the results to be normally distributed).

1) If it is the variance that is of most interest to me then I believe that it is OK to summarise the variance of a group of subjects by taking the mean of their variances?

2) Would it also be meaningful to look at the variance of the subject variances?

3) I believe this (or a similar) statistic will be required if I am to use a T-test to assess whether the difference in variance between groups is statistically significant?

Very simple example (edited to correct the variances as noted in one answer)

Andrew: 5, 2, 7 (var 4.22 6.33)

Bob: 3, 3, 2 (var 0.22 0.33)

Charles: -2, 1 , 0 (var 1.55 2.33)

Diane: 4, 2, 3 (var 0.66 1.00)

Elsa: 6, -4, -1 (var 17.55 26.3)

Fran: 6, 0, 3 (var 6.00 9.00)

Male mean variance (Andrew, Bob and Charles): 2.00 3.00

Female mean variance (Diane, Elsa and Fran): 8.07 12.1

Question: Is there a significant difference between the male and female within-subject variance for this test?

• so can estimate mean and variance of each subject if I understand you well ? Then you make groups of subjects ? Why then not compute mean and variance for each group in the same way as you computed it for a subject (i.e. consider a group as a 'subject') ? Why compare variances with a t-test and not with an F-test ?
– user83346
Sep 13, 2015 at 10:24
• Yes, from the within-subject repeats I can estimate each subject's mean and variance. I then wish to compare the variance for groups of subject (e.g. male vs female) to see if there is a significant difference. I'm not sure I understand your suggestion of considering groups as subjects - can you elaborate? I'm also not sure how an F-test would be preferable to a T-test?
– welf
Sep 14, 2015 at 8:23
• Well, from the within subject repeats you can estimate each subject's mean and variance. I would say that, in a similar way, from the within group repeats you can estimate each group's mean and variance ? If you have the variance for the male and for the female, then you can compare the two variances with an F-test ?
– user83346
Sep 14, 2015 at 8:41
• So for any group, the mean will be the mean of the subject means (equivalent to the mean of all subject tests within the group if all subjects perform the same number of tests?), and the variance will be the variance of the subject means? I feel like this would be disregarding the within-subject variances - hence my suggestion of using the mean of the within-subject variances?
– welf
Sep 14, 2015 at 9:33
• Are you now talking about multiple tests ? because in your question it says "... repeated a test ..." ? I did not say that I would compute the variance of the means ? I would compute the variance of the subjects in the group ?
– user83346
Sep 14, 2015 at 10:31

You can test this by fitting a linear mixed model. A linear mixed model is like a multiple regression model but you can have random effects. The random effects part is needed because you have multiple tests per observation. You will then model score as a function of sex and test, and the subjects are your random effects. Let's enter your test data in R:

subject <- c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6)
score <- c(5,2,7,3,3,2,-2,1,0,4,2,3,6,-4,-1,6,0,3)
test <- rep(c(1,2,3),6)
sex <- c(rep(1,9), rep(0,9))
testdata <- data.frame(subject, score, test, sex)
testdata


The data would look like this:

   subject score test sex
1        1     5    1   1
2        1     2    2   1
3        1     7    3   1
4        2     3    1   1
5        2     3    2   1
6        2     2    3   1
7        3    -2    1   1
8        3     1    2   1
9        3     0    3   1
10       4     4    1   0
11       4     2    2   0
12       4     3    3   0
13       5     6    1   0
14       5    -4    2   0
15       5    -1    3   0
16       6     6    1   0
17       6     0    2   0
18       6     3    3   0


Now we fit two models. Both models use sex and test (1-3 in this case) as fixed effects and subject as random effect. The difference between the models is that in the second model, the variance is allowed to differ between women and men. We then compare the models using the anova() command, and if there is a significant difference, this indicates that the more complex model (the one with the differing variances per sex) provides a better fit, and we thus have indirect evidence that the difference in variance is statistically significant:

m1 <- lme(score ~ sex + factor(test), random=~1|subject, data=testdata)
m2 <- lme(score ~ sex + factor(test), random=~1|subject, weights=varIdent(form=~1|sex), data=testdata)
anova(m1, m2)
Model df      AIC      BIC    logLik   Test    L.Ratio p-value
m1     1  6 167.8734 176.6679 -77.93673
m2     2  7 169.8450 180.1051 -77.92247 1 vs 2 0.02850744  0.8659


There was no difference in this example. But if we change the score for the last female a little (changing score from 3 to -17, increasing the variance) and run the m2 model and the comparison again:

score <- c(5,2,7,3,3,2,-2,1,0,4,2,3,6,-4,-1,6,0,-17)
testdata <- data.frame(subject, score, test, sex)
m2 <- lme(score ~ sex + factor(test), random=~1|subject, weights=varIdent(form=~1|sex), data=testdata)
anova(m1, m2)
Model df      AIC      BIC    logLik   Test  L.Ratio p-value
m1     1  6 105.3948 109.2291 -46.69739
m2     2  7 102.7737 107.2471 -44.38685 1 vs 2 4.621064  0.0316


Now we see a difference in AIC and logLik, and a low p-value which indicates a difference in variance between the sexes.

• This is how it is done. Sep 14, 2015 at 12:38
• IIRC, testing on the edge of the parameter space can provide overly conservative results, and zero variance is definitely at the edge of the parameter space. I'm not sure what the impact would be with the way you've set up your models here (with weights). Sep 14, 2015 at 13:49
• When you use random effects the you implicitly assume a var-covar matrix that is 'compound symmetrie'. The implicit assumption that you make is thus that the correltation between the results of the first and the second test is equal to the one between the second and the third test and the correlation between the second and the third. If that is a reasonable assumption then random effects can be used here, if not you have to use a more general var-covar structure (see my answer)
– user83346
Sep 16, 2015 at 12:00
• Ah, yes of course! I didn't think about the possibility of the subjects becoming better at the test after repeated tries or something like that. I just thought about the tests as a series of unconnected tests (such as questions on an exam or something) that have the same covariance.
– JonB
Sep 16, 2015 at 12:23

The average of variances is probably not the best variable to use.

Variance is a squared result. In many cases, the standard deviation may be the better result to use, because it is on the same scale as your original data. You can then decide on whether use the arithmetic mean, or e.g. root-mean-square again.

Avoid blindly using standard toolbox functions.

The proper variances are:

6.33, 0.33, 2.33, 1, 26.3, 9.


because you must use the unbiased sample variance. Note that these estimates are much higher than yours, because of the small sample size.

As you can see, there is a massive outlier here - 26.3 is way outside your range. Averaging such squared values is not sound. Take the standard deviations instead:

2.52,  0.58,  1.53,  1.00,  5.13,  3.00


The 5.13 is still large, but not as extreme anymore. The mean standard deviation of the males is 1.54, of the females it is 3.04; the average standard deviation of both is $2.29 \pm 1.66$. But you need to be aware that at this sample size, even your estimates of the mean are pretty unreliable.

You have to make informed decision on a number of steps, including:

• variance, or standard deviation? or standard error?
• more than one variance/standard deviation (biased, unbiased)
• more than one mean (arithmetic, harmonic, geometric, power ...)

Any statistic test such as the t-test will come with some assumption on your input data; and doing the wrong choices will have a considerable impact on your result. Sorry, but I cannot save you from studying these differences yourself. The only choice where I'm fairly confident myself is that you need to use the unbiased variant.

For your actual test, you may want to look at the standard error of the mean, too. This may at first appear redundant to the standard deviation, but it is not. It measures how good your estimation of the mean is. It is tighter than the standard deviation; but if you have differences within this standard error, any difference likely is just random.

• Ah, apologies for the use of biassed sample variance, I was trying to put together the example quickly to illustrate my problem. Also thank you for your other useful comments.
– welf
Sep 14, 2015 at 13:03

The fact the same subject does the test multiple times introduces correlation among the observations. The goal is to get an estimate of the var-covar matrix of your data. We will assume that the subjects are independent, but that the score for the same subject on the test are dependent. This means that that var-covar matrix of the data will look like a diagonal matrix, but on the main diagonal you will find $3\times3$ block matrices.

This var-covar matrix can be estimated using the function gls from the package nlme.

Using the same testdata as @Jonas Berge I have the following R-code

library(nlme)
subject <- as.factor(c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6))
score <- c(5,2,7,3,3,2,-2,1,0,4,2,3,6,-4,-1,6,0,3)
test <- (rep(c(1,2,3),6))
sex <- as.factor(c(rep(1,9), rep(0,9)))
testdata <- data.frame(subject, score, test, sex)
testdata

r.1<-gls(score ~ 1+sex,
data=testdata,
corr=corCAR1(form=~ 1|subject),
method="REML")


The corr=corCAR1(form=~ 1|subject) means that I assume an autoregressive pattern in the scores of the test for by subject. You can try with other assumptions on this correlation structure (see help on gls).

Note that we used option method=REML to get unbiased estimates of the variance. (see Why does one have to use REML (instead of ML) for choosing among nested var-covar models?)

With summary(r.1) you see what this gives, but at this point this is not so important. The most important result in this context is that value of the log-likelihood of this model, it will be used as input for a likelihood ratio test further on.

The above model r.1 assumes a correlation structure whereby all subjects have the same variance. In the second stage we will re-estimate a model, but with different variances for the male and female subjects:

r.2<-update(r.1,
weights=varIdent(form=~1|sex),
method="REML")


The option weights=varIdent(form=~1|sex) means that we want different variances by sex. With the summary function one can again find the likelihood for this model.

The second model r.2 has different variances by sex, so it is more general than the first model that assumes that both sexes have the same variance. Therefore the log-likelihood of the second model will be higher than the lig-likelihood of the second model. With a likelihood ratio test (see What are the ''desirable'' statistical properties of the likelihood ratio test?) we can find out whether the difference is significant: the second model r.2 is more general than the first one or the first model is nested within the second one. If we find that the second model's likelihood is significantly higher than the likelihood of the first model, then the model r.2 (with different variances for male and female) fits the data better. So if the likelihood ratio test gives a significant result, then we have reason to believe that the variance differs by sex.

The likelihood ratio test takes two times the difference between the two likelihoods, and this is (asymptotically) a $\chi^2$ with degrees of freedom equal to the difference in number of parameters between teh two models. The second model has one parameter more (two variances, one per sex) so df=1.

The p-value of the test can be found with:

1-pchisq(2*(r.2$logLik - r.1$logLik), df=1)


It is $0.49$ so there is no significant difference in 'fit' between the two models and therefore we keep the simpler model, with male and female having the same variance.

To find the variances one executes summary(r.2). At the bottom you find that the residual standard error is 2.66, squaring this yields the variance, and under 'Variance function' one finds that for sex==1 the variance is 2.66^2 * 1, for sex==0 it is 2.66^2 * 1.61.

• Thank you for your questions and detailed answer. In our tests the subjects are performing a task that is already familiar to them, they get no feedback on their performance in the test, and tests are all conducted within a short space of time; therefore I think it is reasonable to assume that there will not be any correlation in within-subject performance. In this case would you agree with the answer given by @JonB, or would your solution still differ in some way?
– welf
Sep 21, 2015 at 21:10
• @welf: well everyone can make any assumption he likes to make, but what would you do if you estimate $y=\beta_0 + \beta_1x_1 + \beta_2x_2$ then what would you do if you think that $\beta_1$ is zero? Would you estimate $y=\beta_0+\beta_2x_2$ ? Or would you estimate the full equation and perform a hypothesis test to see whether the sample confirms that it could be zero?
– user83346
Sep 22, 2015 at 4:44
• @welf: I would do a hypothesis test, moreover, you might be interested in the value of the two variances, these are easily read from the estimated var-covar matrix while with a mixed effects model you will have to do some ( not so difficult) calculations.
– user83346
Sep 22, 2015 at 4:55
• That makes sense. Out of interest, is there anything explicitly incorrect in my initial suggestion of comparing the mean of within-subject variances? (Or the mean of the standard deviations as @Anony-Mousse mentions)
– welf
Sep 23, 2015 at 21:48

I'm not sure if you can assume normality of the performance variances, but I'm still looking into it. What you're saying makes sense.

Since you are using a within-subjects design (prone to carryover effects) I think it would be most interesting to look at how the variances are changing over time. Are your subjects performing better or worse?

Yes, you would want to use a two-sample Student's t-test to look at that.

• Thank you for your response. As all results per subject are collected over quite a short period of time, in a randomised sequence, we don't expect there to be much change with time. Can you give any more information on your comment about normality?
– welf
Sep 14, 2015 at 8:31

It is not very well known but as an alternative to linear mixed model suggested by Jonas Berge you can perform a F test of equality of variances

https://en.wikipedia.org/wiki/F-test_of_equality_of_variances

• Interesting, I didn't know about that. But if the mean scores (and variance) of the different tests differ between tests, then there is a risk that normality cannot be assumed and it won't be appropriate, according to the wikipedia article. And also, if we just pool the results from all tests for all men vs women, will that not also violate some assumption of independence of observations? So would it the be best to run this f-test once per test?
– JonB
Sep 14, 2015 at 12:41
• This sounds useful, but I'm not sure how to best apply it to my situation, particularly considering the answer given by @Anony-Mousse and the comment above?
– welf
Sep 14, 2015 at 13:49
• well, @Anony-Mousse was right you do need to carefully consider your hypothesis i.e. the thing you want to test. As the F-test of equality of variances aim to test the H0: two groups of variances is equal. Sep 15, 2015 at 10:59
• The F-test can be used to compare the variances, but as I explain in my answer, you have to take into account the dependence structure caused by the repeated measures. With the gls function, you estimate the variance-covariance structure, i.e. also the covariance between repeated measures for the same subject.
– user83346
Sep 16, 2015 at 6:02