Comparing the within-subject variance between two groups of subjects

Question

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?

Answer 1

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.

Answer 2

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.

In your case, I already disagree with your variances.

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.

Answer 3

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.

Answer 4

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.

Answer 5

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