How to interpret verbose output
On 14-05-21 06:15 PM, Stuart Luppescu wrote:
Hello, I'm trying to do a variance decomposition of teacher performance
ratings. Each teacher is observed four times and rated on 9 components.
The model looks like this:
lme8 <- lmer(rating ~
(1|tid.f) + (1|obsorder.f) + (1|comp.f) +
(1|tid.f:comp.f) + (1|obsorder.f:comp.f) + (1|
tid.f:obsorder.f) ,
data=ratings, REML=FALSE, verbose=2)
where tid.f is the teacher identifier, comp.f the component identifier,
and obsorder.f the observation number.
This works fine for the whole dataset, but I want to do it separately by
decile based on the average rating for each teacher, so I added this: ,
subset=ratings$bins==quantile
where the bins are {1, ..., 10} indicating the decile, and quantile is a
loop index used like this: for(quantile in 1:10) {}
It works fine for the lowest decile, but fails for every decile after
that. I get 0.00 for the tid.f variance component, which is the one I'm
really interested in. I have no idea why. I checked the distribution of
average ratings by decile and it all looks unremarkable. I think there
may be a clue in the iteration details as shown by the verbose output
but I don't know how to interpret it. Here it is for decile 1:
This printing comes from minqa::bobyqa. ?bobyqa says:
iprint The value of ?iprint? should be set to ?0, 1, 2 or 3?,
which controls the amount of printing. Specifically, there is
no output if ?iprint=0? and there is output only at the
return if ?iprint=1?. Otherwise, each new value of ?rho? is
printed, with the best vector of variables so far and the
corresponding value of the objective function. Further, each
new value of the objective function with its variables are
output if ?iprint=3?. Default value is ?0?.
npt = 8 , n = 6
Number of points used for quadratic approximation, number of points (parameters)
rhobeg = 0.2 , rhoend = 2e-07 0.020: 16: 18359.8;0.454554 0.735854 0.376302 0.705437 0.787591 0.802885
rho (current trust region radius): number of function evaluations
(?): current value of objective function (PWRSS/deviance); values of
parameters ("theta") -- lower triangle of the Cholesky factor of the RE
variance-covariance matrices, scaled by the residual standard deviation.
0.0020: 32: 18346.2;0.454004 0.699893 0.407403 0.574987 0.702683 0.880645 0.00020: 77: 18263.7;0.432231 0.697895 0.398860 0.0461275 0.326475 1.55039 2.0e-05: 292: 18252.6;0.432051 0.700841 0.394275 0.0472994 0.315551 0.130071 2.0e-06: 339: 18252.6;0.432029 0.700839 0.394515 0.0472285 0.314849 0.129962 2.0e-07: 364: 18252.6;0.432028 0.700822 0.394548 0.0472181 0.314777 0.129904 At return 395: 18252.578: 0.432030 0.700822 0.394554 0.0472204 0.314755 0.129913
print(summary(lme8))
Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: rating ~ (1 | tid.f) + (1 | obsorder.f) + (1 | comp.f) + (1 |
tid.f:comp.f) + (1 | obsorder.f:comp.f) + (1 | tid.f:obsorder.f)
Data: ratings
Subset: ratings$bins == quantile
AIC BIC logLik deviance df.resid
18268.6 18328.1 -9126.3 18252.6 12622
Scaled residuals:
Min 1Q Median 3Q Max
-3.2200 -0.5592 -0.0306 0.5711 4.4942
Random effects:
Groups Name Variance Std.Dev.
tid.f:comp.f (Intercept) 0.033067 0.18184
tid.f:obsorder.f (Intercept) 0.087013 0.29498
tid.f (Intercept) 0.027579 0.16607
obsorder.f:comp.f (Intercept) 0.000395 0.01988
comp.f (Intercept) 0.017551 0.13248
obsorder.f (Intercept) 0.002990 0.05468
Residual 0.177161 0.42090
Number of obs: 12630, groups: tid.f:comp.f, 3159; tid.f:obsorder.f,
1404; tid.f, 351; obsorder.f:comp.f, 45; comp.f, 9; obsorder.f, 5
Fixed effects:
Estimate Std. Error t value
(Intercept) 2.01657 0.05355 37.65
And here for decile 5:
npt = 8 , n = 6
rhobeg = 0.2 , rhoend = 2e-07
0.020: 18: 15847.4;0.261469 0.514852 0.178718 0.0962381
0.152470 0.160530
0.0020: 35: 15749.5;0.356245 0.532056 0.00000 0.0961818
0.219586 0.166497
0.00020: 90: 15742.1;0.349036 0.516678 0.00000 0.0755764
0.330633 0.248861
2.0e-05: 107: 15742.1;0.348612 0.515454 0.00000 0.0766747
0.332641 0.251693
2.0e-06: 250: 15742.1;0.348335 0.515211 0.00000 0.0765898
0.330456 0.260112
2.0e-07: 273: 15742.1;0.348339 0.515207 0.00000 0.0765856
0.330510 0.260162
At return
281: 15742.076: 0.348339 0.515207 1.97296e-07 0.0765858 0.330510
0.260162
print(summary(lme8))
Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: rating ~ (1 | tid.f) + (1 | obsorder.f) + (1 | comp.f) + (1 |
tid.f:comp.f) + (1 | obsorder.f:comp.f) + (1 | tid.f:obsorder.f)
Data: ratings
Subset: ratings$bins == quantile
AIC BIC logLik deviance df.resid
15758.1 15817.7 -7871.0 15742.1 12772
Scaled residuals:
Min 1Q Median 3Q Max
-3.6535 -0.3833 0.1362 0.5225 3.4610
Random effects:
Groups Name Variance Std.Dev.
tid.f:comp.f (Intercept) 0.0192757 0.13884
tid.f:obsorder.f (Intercept) 0.0421666 0.20535
tid.f (Intercept) 0.0000000 0.00000
obsorder.f:comp.f (Intercept) 0.0009318 0.03052
comp.f (Intercept) 0.0173531 0.13173
obsorder.f (Intercept) 0.0107521 0.10369
Residual 0.1588568 0.39857
Number of obs: 12780, groups: tid.f:comp.f, 3195; tid.f:obsorder.f,
1420; tid.f, 355; obsorder.f:comp.f, 45; comp.f, 9; obsorder.f, 5
Fixed effects:
Estimate Std. Error t value
(Intercept) 2.85126 0.06774 42.09
For decile 5, I notice that the column fourth from the right goes to 0
after the first line. Does that mean something? Any ideas about why this
is failing?
Thanks in advance.
It means the second element of the theta vector (which controls the among-'tid.f' variation in the intercept) is going to zero. This is probably the correct numeric value ... I would consider further (1) scaling and centering continuous predictors to see if that helps (2) simulating pseudo-data and fitting it to see how it performs (?simulate.merMod will help with this) (3) using bbmle::slice2D to view 2D slices of the deviance surface.