Should the within-subject variability decrease?

I have a crossover experiment design, detailed as follows.

1. there are 7 sites conducting the same experiment;
2. In the experiment of each sites, 5 different treatments are administered to $n_i, \, (i = 1, 2, \dots, 7)$ subjects during 5 different periods. The $n_i$'s are not necessarily equal;
3. In this experiment, only the fixed treatment effect and random subject effect are of interest. Other effects are ignored.

The variance of the random subject effect is thus the between-subject variability, and the remaining variability from experiment error is within-subject variability.

Three mixed-effects models in terms of different sites are investigated.

(1) Site 1 and 2 combined; (2) Sites other than 1 and 2; (3) All sites.

In each mixed-effect model, fixed treatment effect and random subject effect are specified. After running the three models, I got to know that the within-subject variability for the three are,

(1) 0.7904; (2) 0.8203; (3) 0.8526

respectively.

However, I was told that there results should NOT be like this. Instead, intuitively the within-subject variability from Model (3) should be less than that from Model (2). The argument is that since Model (1) has the smallest within-subject variability, then if we combine the sites (subjects) in Model (1) with those in Model (2), which results to Model (3), we should expect the within-subject variability decrease.

I am not convinced by this the simple intuition. My questions are:

(1) Is the above argument against the results correct?

(2) Are there clearer intuitions that help explain this problem?

(3) Are there any proofs or arguments from theoretical aspect (I mean by derivations) to help understand this?

EDIT:

My description to the problem might not be quite clear. But you can think of subjects from site 1 and 2 are of high quality control since the Model (1) produces lower within-subject variability. Thus the rational which is against the current numerical results is, if these "better" subjects are combined with others (Site 3 to 7), it's expected the within-subject variability decrease. Thus the within-subject variability from Model (3) should be smaller than that of Model (2). What I am asking is whether this rational necessarily holds. If so, how? If not, then why?

END EDIT

EDIT-2:

I am posting some "made-up" data for anyone who might be curious about. It's of longitudinal type. There are four variables "Subject", "Treatment", "Site", and response variable "y". As I said above, each site has some subjects for the experiment. There are 5 different treatment levels, denoted as A, B, C, D, and E.

The model is pretty simple, as shown in following sample SAS code, with fixed treatment effect and random subject effect. Intentionally, I am not considering other effects (like site effect). The code applies to all three models, except that the "data = case2" corresponds to different (sub)sets of data with different sites included.

proc mixed data = case2 asycov cl covtest plots;
class Subject   Treatment;
model y= Treatment/solution;
random Subject;
lsmeans TRT/pdiff cl;
run;

Here is the sample data.

Subject Treatment   Site    y
1-1     D              1    5.68387
1-1     E              1    5.65
1-1     C              1    4.45098
1-1     A              1    0.79048
1-1     B              1    4.50455
1-2     C              1    4.13208
1-2     D              1    5.10459
1-2     B              1    4.34468
1-2     E              1    5.07556
1-2     A              1    2.36296
1-3     B              1    -0.77037
1-3     C              1    0.59167
1-3     A              1    -1.53191
1-3     D              1    3.42
1-3     E              1    2.89231
...
2-1     D              2    4.80312
2-1     E              2    5.60606
2-1     C              2    5.38
2-1     A              2    0.39474
2-1     B              2    3.97714
2-2     C              2    4.46667
2-2     D              2    5.73333
2-2     B              2    6.17391
2-2     E              2    4.86957
2-2     A              2    -0.02
...
...
7-1     E              7    4.4875
7-1     A              7    -2.10667
7-1     D              7    5.71579
7-1     B              7    1.17895
7-1     C              7    0.97576
7-2     D              7    3.05083
7-2     E              7    3.25473
7-2     C              7    2.7925
7-2     A              7    1.23
7-2     B              7    3.98769
7-3     C              7    4.32754
7-3     D              7    5.2875
7-3     B              7    4.46575
7-3     E              7    4.12787
7-3     A              7    0.61481

END EDIT-2

• Do the models include site effects (either fixed or random)? Or are sites ignored (maybe that is what you mean by "combined")? – Jake Westfall Nov 13 '14 at 7:20
• @JakeWestfall, yes, there is no site effect in all the three models. Only fixed treatment effect and random subject effect are considered. – Aaron Zeng Nov 13 '14 at 13:18
• @AaronZeng Perhaps I misunderstand the situation... but I would think that if you are trying to fit the same model to more data the model fit to each individual subject is likely to be worse. Thus more variance will be attributed to "within-subjects" which acts as a "variance sink" for the model. – Livid Nov 13 '14 at 20:07
• @Livid, I am sure what you mean by "fit to each individual subject". But here is my situation. You can think of subjects from site 1 and 2 are of high quality control due to lower within-subject variability. Thus the rational is, if these "better" subjects are combined with others (Site 3 to 7), it's expected the within-subject variability decrease. I am asking whether this rational necessarily holds. If so, how? If not, then why? – Aaron Zeng Nov 13 '14 at 20:22
• @AaronZeng I think I know what is going on here, if I am correct then within-subject variance will increase with the number of sites (what I called "subjects"). Can you provide a table of some example data (it can be "made up" but similar to the real data) and code (R if possible) for the analysis you are performing? – Livid Nov 13 '14 at 21:29

The situation is somewhat different from what I thought. In case someone else can answer it, here is the data and my understanding of (I do not know SAS well...) the analysis in R:

require(beanplot)
require(lme4)

dat.full<-structure(list(Subject = c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L,
6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L),
Treatment = structure(c(4L, 5L, 3L, 1L, 2L, 3L, 4L, 2L, 5L,
1L, 2L, 3L, 1L, 4L, 5L, 4L, 5L, 3L, 1L, 2L, 3L, 4L, 2L, 5L,
1L, 5L, 1L, 4L, 2L, 3L, 4L, 5L, 3L, 1L, 2L, 3L, 4L, 2L, 5L,
1L), .Label = c("A", "B", "C", "D", "E"), class = "factor"),
Site = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 7L, 7L,
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L), y = c(5.68387,
5.65, 4.45098, 0.79048, 4.50455, 4.13208, 5.10459, 4.34468,
5.07556, 2.36296, -0.77037, 0.59167, -1.53191, 3.42, 2.89231,
4.80312, 5.60606, 5.38, 0.39474, 3.97714, 4.46667, 5.73333,
6.17391, 4.86957, -0.02, 4.4875, -2.10667, 5.71579, 1.17895,
0.97576, 3.05083, 3.25473, 2.7925, 1.23, 3.98769, 4.32754,
5.2875, 4.46575, 4.12787, 0.61481)), .Names = c("Subject",
"Treatment", "Site", "y"), class = "data.frame", row.names = c(NA,
-40L))

par(mfrow=c(2,2))
beanplot(dat.full$y~dat.full$Subject, xlab="Subject",
ylim=c(min(dat.full$y)-1, max(dat.full$y)+1))
beanplot(dat.full$y~dat.full$Site, xlab="Site",
ylim=c(min(dat.full$y)-1, max(dat.full$y)+1))
beanplot(dat.full$y~dat.full$Treatment, xlab="Treatment",
ylim=c(min(dat.full$y)-1, max(dat.full$y)+1))

Nsubj<-length(unique(dat.full$Subject)) plot(0,0, type="n",xlim=c(1,5), xaxt="n", xlab="Treatment", ylab="Y", ylim=c(min(dat.full$y)-1, max(dat.full$y)+1)) axis(side=1, at=as.numeric(unique(dat.full$Treatment)),
labels=unique(dat.full$Treatment)) for(i in 1:Nsubj){ sub<-dat.full[which(dat.full$Subject==unique(dat.full$Subject)[i]),] sub<-sub[sort(as.numeric(sub$Treatment), index.return=T)$ix,] lines(sub$Treatment,sub$y, col=rainbow(Nsubj)[i], lwd=3) } legend("bottomright", legend=unique(dat.full$Subject), col=rainbow(Nsubj),
lwd=4, title="Subject", ncol=2)

dat1<-dat.full[which(dat.full$Site %in% 1),] dat27<-dat.full[which(dat.full$Site %in% c(2,7)),]

fit.all<-lmer(dat.full$y~dat.full$Treatment + (1|dat.full$Subject)) fit1<-lmer(dat1$y~dat1$Treatment + (1|dat1$Subject))
fit27<-lmer(dat27$y~dat27$Treatment + (1|dat27$Subject)) > fit.all Linear mixed model fit by REML ['lmerMod'] Formula: dat.full$y ~ dat.full$Treatment + (1 | dat.full$Subject)
REML criterion at convergence: 127.467
Random effects:
Groups           Name        Std.Dev.
dat.full$Subject (Intercept) 1.142 Residual 1.064 Number of obs: 40, groups: dat.full$Subject, 8
Fixed Effects:
(Intercept)  dat.full$TreatmentB dat.full$TreatmentC
0.2168               3.2660               3.1728
dat.full$TreatmentD dat.full$TreatmentE
4.6331               4.2786
> fit1
Linear mixed model fit by REML ['lmerMod']
Formula: dat1$y ~ dat1$Treatment + (1 | dat1$Subject) REML criterion at convergence: 37.0667 Random effects: Groups Name Std.Dev. dat1$Subject (Intercept) 1.8608
Residual                 0.8506
Number of obs: 15, groups:  dat1$Subject, 3 Fixed Effects: (Intercept) dat1$TreatmentB  dat1$TreatmentC dat1$TreatmentD
0.5405           2.1524           2.5177           4.1956
dat1$TreatmentE 3.9988 > fit27 Linear mixed model fit by REML ['lmerMod'] Formula: dat27$y ~ dat27$Treatment + (1 | dat27$Subject)
REML criterion at convergence: 76.0409
Random effects:
Groups        Name        Std.Dev.
dat27$Subject (Intercept) 0.7462 Residual 1.1876 Number of obs: 25, groups: dat27$Subject, 5
Fixed Effects:
(Intercept)  dat27$TreatmentB dat27$TreatmentC  dat27$TreatmentD 0.02258 3.93411 3.56592 4.89554 dat27$TreatmentE
4.44657

My understanding is that @AaronZeng wants to know under what conditions the residuals can be fit.all>fit27>fit1. This is not the case in the example data.