Hand-calculation of predicted values with splines

I'm trying to hand-calculate the predicted values from a simple model
that includes a spline. To illustrate, I'll the orange tree data for
tree 1.

tree1 <- subset(orange, Tree==1)
m1 <- asreml(circ ~ 1 + lin(x) , random= ~ spl(x), data=tree1)
m1 <- update(m1)
p1 <- predict(m1, classify="x",
              predictpoints=list(x=seq(from=118, to=1582, length=20)))
p1 <- p1$predictions$x$pvals
# Plot data and predicted value
plot(circ~x, data=tree1, type="l")
text(1000,40,"Observed", col="black")
with(p1, points(predicted.value~x, type="l", col="blue"))
text(1000,45, "Spline predicted", col="blue")
# Plot the fixed linear part of the prediction
avgx <- mean(tree1$x)
intx <- coef(m1)$fixed[2]
slopex <- coef(m1)$fixed[1]
points(tree1$x, intx + slopex * tree1$x, col="red", type='b')
text(1000, 50, "Linear predicted", col="red")

An ordinary (non mixed model) cubic spline is of the form
 Si(x) = ai (xi - x)^3 + bi(xi - x)^2 + ci(xi - x) + di
but for this mixed model maybe we need to sweep out (subtract) the linear part?

# Observed - linear
yvals <- tree1$circ - (intx + slopex * tree1$x)

# Now try to calculate the predicted value at 888.5263
sp <- function(newx, xi, yi, coefs){
  hi <- diff(xi)
  di <- yvals
  bi <- c(0, coefs, 0) # b0 and bn are fixed at 0
  ai <- diff(bi)/(3*hi)
  ci <- diff(yi)/hi - (bi[-1] + 2*bi[-length(bi)])/(3*hi)
  # In which interval is the new point?
  i <- max(which(newx >= xi))
  ai[i]*(newx-xi[3])^3 + bi[i]*(newx-xi[3])^2 + ci[i]*(newx-xi[3]) + di[i]
}
# Intercept + slope*888.5263 + spline at 888.5263
intx + slopex*888.5263 + sp(888.5263, tree1$x, yvals, coef(m1)$random)

This gives a value of -55642, so obviously something is wrong. I'm
not sure that I'm interpreting the "fixed linear, random spline" idea
correctly, I'm not sure that I'm using the spline coefficients
correctly and I'm not sure that I have parameterized the spline in the
same way that asreml does.

Any clarification (or other worked examples) would be appreciated.

Kevin Wright
Received on Sat Nov 06 2007 - 17:01:44 EST