Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

I'm trying to analyze effect of Year on variable logInd for particular group of individuals (I have 3 groups). The simplest model:

> fix1 = lm(logInd ~ 0 + Group + Year:Group, data = mydata)
> summary(fix1)

Call:
lm(formula = logInd ~ 0 + Group + Year:Group, data = mydata)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.5835 -0.3543 -0.0024  0.3944  4.7294 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
Group1       4.6395740  0.0466217  99.515  < 2e-16 ***
Group2       4.8094268  0.0534118  90.044  < 2e-16 ***
Group3       4.5607287  0.0561066  81.287  < 2e-16 ***
Group1:Year -0.0084165  0.0027144  -3.101  0.00195 ** 
Group2:Year  0.0032369  0.0031098   1.041  0.29802    
Group3:Year  0.0006081  0.0032666   0.186  0.85235    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.7926 on 2981 degrees of freedom
Multiple R-squared: 0.9717,     Adjusted R-squared: 0.9716 
F-statistic: 1.705e+04 on 6 and 2981 DF,  p-value: < 2.2e-16 

We can see the Group1 is significantly declining, the Groups2 and 3 increasing but not significantly so.

Clearly the individual should be random effect, so I introduce random intercept effect for each individual:

> mix1a = lmer(logInd ~ 0 + Group + Year:Group + (1|Individual), data = mydata)
> summary(mix1a)
Linear mixed model fit by REML 
Formula: logInd ~ 0 + Group + Year:Group + (1 | Individual) 
   Data: mydata 
  AIC  BIC logLik deviance REMLdev
 4727 4775  -2356     4671    4711
Random effects:
 Groups     Name        Variance Std.Dev.
 Individual (Intercept) 0.39357  0.62735 
 Residual               0.24532  0.49530 
Number of obs: 2987, groups: Individual, 103

Fixed effects:
              Estimate Std. Error t value
Group1       4.6395740  0.1010868   45.90
Group2       4.8094268  0.1158095   41.53
Group3       4.5607287  0.1216522   37.49
Group1:Year -0.0084165  0.0016963   -4.96
Group2:Year  0.0032369  0.0019433    1.67
Group3:Year  0.0006081  0.0020414    0.30

Correlation of Fixed Effects:
            Group1 Group2 Group3 Grp1:Y Grp2:Y
Group2       0.000                            
Group3       0.000  0.000                     
Group1:Year -0.252  0.000  0.000              
Group2:Year  0.000 -0.252  0.000  0.000       
Group3:Year  0.000  0.000 -0.252  0.000  0.000

It had an expected effect - the SE of slopes (coefficients Group1-3:Year) are now lower and the residual SE is also lower.

The individuals are also different in slope so I also introduced the random slope effect:

> mix1c = lmer(logInd ~ 0 + Group + Year:Group + (1 + Year|Individual), data = mydata)
> summary(mix1c)
Linear mixed model fit by REML 
Formula: logInd ~ 0 + Group + Year:Group + (1 + Year | Individual) 
   Data: mydata 
  AIC  BIC logLik deviance REMLdev
 2941 3001  -1461     2885    2921
Random effects:
 Groups     Name        Variance  Std.Dev. Corr   
 Individual (Intercept) 0.1054790 0.324775        
            Year        0.0017447 0.041769 -0.246 
 Residual               0.1223920 0.349846        
Number of obs: 2987, groups: Individual, 103

Fixed effects:
              Estimate Std. Error t value
Group1       4.6395740  0.0541746   85.64
Group2       4.8094268  0.0620648   77.49
Group3       4.5607287  0.0651960   69.95
Group1:Year -0.0084165  0.0065557   -1.28
Group2:Year  0.0032369  0.0075105    0.43
Group3:Year  0.0006081  0.0078894    0.08

Correlation of Fixed Effects:
            Group1 Group2 Group3 Grp1:Y Grp2:Y
Group2       0.000                            
Group3       0.000  0.000                     
Group1:Year -0.285  0.000  0.000              
Group2:Year  0.000 -0.285  0.000  0.000       
Group3:Year  0.000  0.000 -0.285  0.000  0.000

But now, contrary to the expectation, the SE of slopes (coefficients Group1-3:Year) are now much higher, even higher than with no random effect at all!

How is this possible? I would expect that the random effect will "eat" the unexplained variability and increase "sureness" of the estimate!

However, the residual SE behaves as expected - it is lower than in the random intercept model.

Here is the data if needed.

Edit

Now I realized astonishing fact. If I do the linear regression for each individual separately and then run ANOVA on the resultant slopes, I get exactly the same result as the random slope model! Would you know why?

indivSlope = c()
for (indiv in 1:103) {
    mod1 = lm(logInd ~ Year, data = mydata[mydata$Individual == indiv,])
    indivSlope[indiv] = coef(mod1)['Year']
}

indivGroup = unique(mydata[,c("Individual", "Group")])[,"Group"]


anova1 = lm(indivSlope ~ 0 + indivGroup)
summary(anova1)

Call:
lm(formula = indivSlope ~ 0 + indivGroup)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.176288 -0.016502  0.004692  0.020316  0.153086 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
indivGroup1 -0.0084165  0.0065555  -1.284    0.202
indivGroup2  0.0032369  0.0075103   0.431    0.667
indivGroup3  0.0006081  0.0078892   0.077    0.939

Residual standard error: 0.04248 on 100 degrees of freedom
Multiple R-squared: 0.01807,    Adjusted R-squared: -0.01139 
F-statistic: 0.6133 on 3 and 100 DF,  p-value: 0.6079 

Here is the data if needed.

share|improve this question
    
You need a year fixed effect if you're going to have a year:group interaction fixed effect. In general, you can't include an interaction term without also including the main effects. Do you really think there's no fixed component to the year effect? And, if so, how could there be a fixed year:group interaction? –  John Jul 10 '12 at 19:57
    
And, why no fixed intercept? You can have both, fixed and random. –  John Jul 10 '12 at 19:59
    
@John, this model is completely valid. This is only an issue of the desired coding of the categorical variable. This way Group$i$ is the intercept of Group $i$, and Group$i$:Year is the slope within Group $i$. If the main effect of Year and the intercept are included, then the estimates would be the differences of the intercept of Group $i$ and Group 1, and similarly with slopes. –  Aniko Jul 10 '12 at 20:02
    
@John, this is off-topic to my question, nevertheless: believe me, this is OK, I did a lot experiments with that. My first lm model is entirely equivalent to logInd ~ Year*Group, only the coefficients are in different shape, nothing more. Depends on your taste and what shape of coefficients you like, nothing more. There's no exclusion of "Year main effect" in my 1st model as you write... logInd ~ Year*Group does exactly the same, the Year coefficient then is not the main effect, but the Group1:Year. –  Curious Jul 10 '12 at 20:04
    
OK, neat, hadn't considered both the 0 intercept and Group being categorical. –  John Jul 10 '12 at 20:20
add comment

1 Answer 1

up vote 8 down vote accepted

I think the problem is with your expectations:) Note that when you added a random intercept for each individual, the standard error of the intercepts increased. Since each individual can have his/her own intercept, the group average is less certain. The same thing happened with the random slope: you are not estimating one common (within-group) slope anymore, but the average of varying slopes.

EDIT: Why doesn't a better model give a more precise estimate?

Let's think about it the other way around: why does the initial model underestimate the standard error? It assumes independence of observations that are not independent. The second model relaxes that assumption (in a way that affects the intercepts), and the third relaxes it further.

EDIT 2: relationship with many patient-specific models

Your observation is a known property (and if you had only two years, then the random effects model would be equivalent to a paired t-test). I don't think I can manage a real proof, but perhaps writing out the two models will make the relationship clearer. Let's ignore the grouping variable, as it would just complicate the notation. I will use greek letters for random effects, and latin letters for fixed effects.

The random effects model is ($i$ - subject, $j$ - replicate within subject): $$Y_{ij} = a + \alpha_i + (b+\beta_i)x_{ij} + \epsilon_{ij},$$ where $(\alpha_i,\beta_i)'\sim N(0,\Sigma)$ and $\epsilon_{ij}\sim N(0,\sigma^2)$.

When you fit separate models for each subject, then $$Y_{ij} = a_i + b_i x_{ij}+ \epsilon_{ij},$$ where $\epsilon_{ij}\sim N(0,\sigma_i^2)$.

[Note: the following is really just handwaving:]

You can see a lot of similarities between these two models with $a_i$ corresponding to $a+\alpha_i$ and $b_i$ to $b+\beta_i$. The average of $b_i$'s corresponds to $b$, because the random effects average to 0. The unconstrained correlation of the random intercept and slope leads to the fact that the models just can be fitted separately. I am not sure how the single $\sigma$ assumption meshes with the subject-specific $\sigma_i$, but I would assume that $\alpha_i$ picks up the difference.

share|improve this answer
    
Thanks Aniko. You are right, my computations confirm that, but I would like to see why... It seems contra-intuitive. I improved the models - by introducing random effects I described the error structure better. Residual error confirms it - is lower and lower. So with these better, more precise models I would expect more precise slope... I know I'm wrong somewhere, please help me see it. –  Curious Jul 10 '12 at 19:30
    
Thanks Aniko, that's an interesting point of view! I am only interested in slopes (Group*:Year), not intercept here.. so my first step of introducing random itcept effect relaxed that independence assumption and lead to the lower SE.. (of slope..) and then the next step was probably too much (??) and did the contrary (even worse SE..) .. maybe I need to think about it, thanks. –  Curious Jul 10 '12 at 20:15
    
Now I'm also astonished by very interesting fact - please see my edit. Would you know why is that? –  Curious Jul 10 '12 at 20:15
    
I don't think the independence assumption was relaxed too much! It was wrong to begin with. –  Aniko Jul 10 '12 at 20:16
3  
Tomas, a "precise" model does not mean the estimates will be more precise. As an extreme example, take any data-free model you like, such as one that predicts all responses are zero. This model is absolutely certain in its estimate of zero. It is therefore as precise as one could possibly get--but it's probably as wrong as possible, too. Giving a model greater scope to fit parameters therefore usually means those parameters are fit with less precision, not more. A better model, because it can quantify the uncertainty not captured by a worse model, often has larger standard errors. –  whuber Jul 10 '12 at 20:17
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.