From: Kevin Wright <kw.statr_at_GMAIL.COM>

Date: Tue, 6 Nov 2007 17:01:44 -0600

Date: Tue, 6 Nov 2007 17:01:44 -0600

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

This webpage is part of the ASReml-l discussion list archives 2004-2010. More information on ASReml can be found at the VSN website. This discussion list is now deprecated - please use the VSN forum for discussion on ASReml. (These online archives were generated using the hypermail package.)