[R-sig-ME] lme random effects in additive models with interaction

[elifnurdogruoz at mynet.com]

Thanks for the warning, I hope it seems better now.

My data is like following and there is a response variable y:

Time    Size            Charge          Density   Replication
3       small             +               low      1
.        .                .                .       .
.        .                .                .       .
9       small             +               low      1
3       big               +               low      .
.        .                .                .       .
.        .                .                .       .
9       big               +               low      1
3       small             -               low      1
.        .                .                .       .
.        .                .                .       .
9       small             -               low      1
3       big               -               low      1
.        .                .                .       .
.        .                .                .       .
9       big               -               low      1
3       small             +               high     1
.        .                .                .       .
.        .                .                .       .
9       small             +               high     1
3       big               +               high     1
.        .                .                .       .
.        .                .                .       .
9       big               +               high     1
3       small             -               high     1
.        .                .                .       .
.        .                .                .       .
9       small             -               high     1
3       big               -               high     1
.        .                .                .       .
.        .                .                .       .
9       big               -               high     1
3       small             +               low      2
.        .                .                .       .
.        .                .                .       .
9       small             +               low      2
3       big               +               low      2
.        .                .                .       .
.        .                .                .       .
9       big               +               low      2
3       small             -               low      2
.        .                .                .       .
.        .                .                .       .
9       small             -               low      2
3       big               -               low      2
.        .                .                .       .
.        .                .                .       .
9       big               -               low      2
3       small             +               high     2
.        .                .                .       .
.        .                .                .       .
9       small             +               high     2
3       big               +               high     2
.        .                .                .       .
.        .                .                .       .
9       big               +               high     2
3       small             -               high     2
.        .                .                .       .
.        .                .                .       .
9       small             -               high     2
3       big               -               high     2
.        .                .                .       .
.        .                .                .       .
9       big               -               high     2


My code with comments:

##this function selects the knots
default.knots <- function(x,num.knots)
{
if (missing(num.knots))
num.knots <- max(5,min(floor(length(unique(x))/4),35))
return(quantile(unique(x),seq(0,1,length=
(num.knots+2))[-c(1,(num.knots+2))]))
}

knots <- default.knots(Time)

z <- outer(Time, knots, "-")
z <- z * (z > 0)
z<-z^2

i.size50 <- I(Size==small)
i.chargepos <- I(Charge=="+")
i.densitylow <- I(Density==low)

##Create X and Z matrices,  I put interactions because I want
intercept to be zero at time 0.
X <- cbind( I(Time^2),Time*i.size50,Time*i.chargepos,Time*i.densitylow)
Z <- cbind( z, z*i.size50, z*i.chargepos,z*i.densitylow)


K <- length(knots)

## form blocked diagonal matrix Z to specify which columns of Z are
used for each group

block.ind <- list(1:K, (K+1):(2*K),(2*K+1):(3*K),(3*K+1):(4*K))
Z.block <- list()
for (i in 1:length(block.ind))
Z.block[[i]] <-
as.formula(paste("~Z[,c(",paste(block.ind[[i]],collapse=","),")]-1"))

##create dummy grouping variable since groupedData object is required for lme
group <- rep(1, length(Time))
model.data <- groupedData(y~X|group, data=data.frame(X, y))

fit <- lme(y~-1+X, data=model.data, random=pdBlocked(list(
pdBlocked(Z.block,pdClass="pdIdent"), pdIdent(~-1+ Replication) ) ) )

It gives error : "Error:  getResponseFormula(el) : "Form" must be a
two sided formula"

I read that thread and I know that 2 replication is not enough to be
modelled as random effect.
However, I want to have the correct model even if I get imprecise
estimation of variance.
I know I have a mistake in forming positive definite matrices in random
part, but I couldn't figure out what is wrong.
Thanks in advance,