> library(MASS)
> data(petrol)
> attach(petrol)
>
> library(nlme)
Loading required package: lattice
Loading required package: grid
> out <- lme(Y ~ SG + VP + V10 + EP, random = ~ 1 | No, method = "ML")
> out.middle <- lme(Y ~ V10 + EP, random = ~ 1 | No, method = "ML")
> out.little <- lme(Y ~ EP, random = ~ 1 | No, method = "ML")
> anova(out.little, out.middle, out)
           Model df      AIC      BIC    logLik   Test  L.Ratio p-value
out.little     1  4 178.6625 184.5255 -85.33126
out.middle     2  5 149.6119 156.9406 -69.80594 1 vs 2 31.05064  <.0001
out            3  7 149.3833 159.6435 -67.69166 2 vs 3  4.22855  0.1207
>  ### lets go with out.middle.  Go back to REML
>  out <- lme(Y ~ V10 + EP, random = ~ 1 | No)
>  summary(out)
Linear mixed-effects model fit by REML
 Data: NULL
       AIC     BIC    logLik
  163.9555 170.792 -76.97775

Random effects:
 Formula: ~1 | No
        (Intercept) Residual
StdDev:    1.607878 1.868132

Fixed effects: Y ~ V10 + EP
                Value Std.Error DF    t-value p-value
(Intercept) 18.209782  4.109272 21   4.431389   2e-04
V10         -0.210971  0.017005  8 -12.406393  <.0001
EP           0.157757  0.005555 21  28.396757  <.0001
 Correlation:
    (Intr) V10
V10 -0.888
EP  -0.185 -0.265

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max
-1.7580007 -0.4209443 -0.1503591  0.4460961  1.7842546

Number of Observations: 32
Number of Groups: 10
> ### BLUE
>  fixed.effects(out)
(Intercept)         V10          EP
 18.2097822  -0.2109710   0.1577568
>  ### BLUP
>  coef(out)
  (Intercept)       V10        EP
A    19.71307 -0.210971 0.1577568
B    16.78038 -0.210971 0.1577568
C    20.21469 -0.210971 0.1577568
D    15.93946 -0.210971 0.1577568
E    17.98715 -0.210971 0.1577568
F    17.95837 -0.210971 0.1577568
G    18.77394 -0.210971 0.1577568
H    18.41735 -0.210971 0.1577568
I    17.71737 -0.210971 0.1577568
J    18.59605 -0.210971 0.1577568
>
>  foo <- data.frame(No = petrol$No[1], SG = 35, VP = 3.1, V10 = 250, EP = 300)
> foo
  No SG  VP V10  EP
1  A 35 3.1 250 300
> predict(out, newdata = foo, level = 0)
[1] 12.79407
attr(,"label")
[1] "Predicted values"
>  predict(out, newdata = foo, level = 1)
Error in val[revOrder, level + 1] : incorrect number of dimensions
>