# Why does fixed effect p-value in a mixed model act unintuitively?

I have a problem with calculating significance of a fixed effect in simple nested design experiment. Suppose we have a simple dataset:

individual  class   ind_mean    value   replicate
1   A   -0.37651119 -1.23495005 1
1   A   -0.37651119 1.24107015  2
1   A   -0.37651119 -0.41977074 3
1   A   -0.37651119 -1.09239410 4
2   A   0.01519211  0.15858308  1
2   A   0.01519211  0.69513257  2
2   A   0.01519211  0.05390856  3
2   A   0.01519211  -0.84685576 4
3   A   -0.19937009 0.49203770  1
3   A   -0.19937009 0.06218100  2
3   A   -0.19937009 -1.12546256 3
3   A   -0.19937009 -0.22623652 4
4   A   0.62853718  1.63792462  1
4   A   0.62853718  1.26212834  2
4   A   0.62853718  -0.52029892 3
4   A   0.62853718  0.13439466  4
5   B   2.38620939  2.67706019  1
5   B   2.38620939  1.71849546  2
5   B   2.38620939  2.31713740  3
5   B   2.38620939  2.83214451  4
6   B   0.91834030  -0.13221319 1
6   B   0.91834030  2.41841628  2
6   B   0.91834030  0.25367274  3
6   B   0.91834030  1.13348538  4
7   B   2.12694664  2.22003887  1
7   B   2.12694664  1.31492118  2
7   B   2.12694664  2.53495794  3
7   B   2.12694664  2.43786856  4
8   B   1.52867659  0.65959322  1
8   B   1.52867659  1.23188018  2
8   B   1.52867659  1.58976692  3
8   B   1.52867659  2.63346604  4


Now if I fit the lme model to this data, I get

> summary(lme(value ~ class, .~1|individual, data = d))
...
Random effects:
Formula: . ~ 1 | individual
(Intercept)  Residual
StdDev:   0.3643563 0.8430941

Fixed effects: value ~ class
Value Std.Error DF  t-value p-value
(Intercept) 0.016962 0.2785935 24 0.060884  0.9520
classB      1.723081 0.3939908  6 4.373405  0.0047
...


Everything seems to be fine. However, I get exactly the same fixed parameter estimates, statistics and p-values if I would fit the model on the column ind_mean that instead of the original values has the averages for each individual.

> summary(lme(ind_mean ~ class, .~1|individual, data = d))
...
Random effects:
Formula: . ~ 1 | individual
(Intercept)     Residual
StdDev:   0.5571871 1.854137e-16

Fixed effects: ind_mean ~ class
Value Std.Error DF  t-value p-value
(Intercept) 0.016962 0.2785935 24 0.060884  0.9520
classB      1.723081 0.3939908  6 4.373405  0.0047
...


Actually this p-value coincides with the ordinary t-test p-value, where I would run the comparison on the averages of individuals. This of course makes sense, since the additional observations in this case do not add any information. However, in the case of original model the additional observations do add information. In fact, the data was generated without any individual influence at all, so the correct statistics would be.

> summary(lm(value ~ class, data = d))
...
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.01696    0.22597   0.075    0.941
classB       1.72308    0.31957   5.392  7.7e-06 ***
...


My statistical intuition would suggest that the p-value from the first model should be somewhere between the ones of third and second model, not exactly the same as in the second model. Is my intuition wrong? Am I using the method in a wrong way? Are there some methods that would give me more appropriate p-values?

$$Y_{ijk}=\beta_0+\beta_i+u_{ij}+e_{ijk}$$
Using the properties of the errors, $e_{ijk}$, The corresponding model for the means is given as:
$$\overline{Y}_{ij}=\beta_0+\beta_i+u_{ij}+\overline{e}_{ij}$$
where the variance of $\overline{e}_{ij}$ is now $\frac{\sigma^2_e}{n_{ij}}$, where $n_{ij}=4$ in your data set (ie "balanced"). Now this shows that the individual means are sufficient for the betas, given the variance parameters. This is why they dont matter. It happens with any divisible distribution that has non unique predictors.
Note that the variance components are $0.55=\sqrt{\sigma^2_u+\frac{\sigma^2_e}{4}}$, which is the variance of $u_{ij}+\overline{e}_{ij}$.