## Required packages
require(ggplot2)
require(lme4)
require(lcmm)
require(gam)

## Simulate and plot raw data
npts <- 50
ntime <- 10
dt <- data.frame(x=rep(seq(0,1,len=ntime), npts),
                 id=as.factor(sort(rep(1:npts, ntime))))
dt$randint <- rnorm(npts, mean=2, sd=0.02)[dt$id]
dt$mu <- dt$randint + 0.7*dt$x + 2*dt$x^2 - 2.4*dt$x^3
dt$y <- rnorm(nrow(dt), mean=dt$mu, sd=0.03)
dt$id <- as.factor(dt$id)
(pp <- ggplot(data=dt, aes(x, y, group=id, color=id))
    + geom_line(col=grey(0.7), show.legend=TRUE))

## Fit spline with lme4 (comparator)
lm1 <- lmer(y~ns(x, df=5) + (1|id), dat=dt, REML=FALSE)

## Fit spline with hlme (ng=1 for simplicity)
hm1 <- hlme(y~ns(x, df=5),subject="id",ng=1,data=dt)

## Compare predictions from both models when using a subset
## of the x used to fit the models (Here is the problem)
newdt <- dt[1:ntime, c("id", "x")]
newdt2<- dt[1:5, c("id", "x")]
newdt2$plm1 <- predict(lm1, newdata=newdt2, re.form=~0)
newdt2$phm1 <- as.numeric(predictY(hm1, newdt2, var.time="x")$pred)
newdt2$smeans <- tapply(dt$y, dt$x, mean)[1:5]
(pp + geom_line(data=newdt2, aes(x=x, y=plm1, color="a"), size=2, show.legend=TRUE)
    + geom_line(data=newdt2, aes(x=x, y=phm1, color="b"), linetype=2, size=2, show.legend=TRUE)
    + geom_point(data=newdt2, aes(x=x, y=smeans, color="c"), size=2, show.legend=TRUE)
    + ylab("y")
    + scale_color_manual(name = "",
                         values = c("a" = "darkblue", "b" = "red", "c"="yellow"),
                         labels=c("predict.merMod(lme4)","predictY(lcmm)","Sample means"))
    + guides(colour = guide_legend(override.aes = list(linetype = c(1,2,0), shape=c(NA,NA,19)))))


## 1) I tried to explicitly provide knots and Boundary.knots but get an error
kk <- quantile(dt$x, prob=seq(0,1,len=5))[2:4]
rr <- range(dt$x)
hm1 <- hlme(y~ns(x,knots=kk, Boundary.knots=rr),subject="id",ng=1,data=dt)

## but get the following Error:
## Be patient, hlme is running ...
## Error in (function (..., row.names = NULL, check.rows = FALSE, check.names = TRUE,  :
##   arguments imply differing number of rows: 500, 3, 2


## 2) Thus tried  to explicitly provide the basis but get an error
mybasis<- ns(dt$x,knots=kk, Boundary.knots=rr)
hm1 <- hlme(y~mybasis,subject="id",ng=1,data=dt)
mynewbasis <- predict(mybasis, newdt2$x)
predictY(hm1, newdata=mynewbasis, var.time="x")

## but get the following Error:
## newdata should at least include the following covariates:
## mybasis NA
## Error in predictY.hlme(hm1, newdata = mynewbasis, var.time = "x") :
##   see above