linear mixed-effects models with two random variables?

Dear R users,
Does anyone knows how to run a glmm with one fixed factor and 2 random
numeric variables (indices)? Is there any way to force in the model a
separate interaction of those random variables with the fixed one?
I hope you can help me.

#eg.
Reserve <- rep(c("In","Out"), 100)
fReserve <- factor(Reserve)
DivBoulders <- rep (c(1.23,2.4,1.26,1.78,1.97,1.35,1.23,2.4,1.26,1.78), 20)
Roughness <- rep(c(3.45,2.56,1.32,5.67,3.73,3.57,2.66,1.52,7.67,2.73),20)
Biomass <- rep(c(8,5.3,3.5,12,25.4,10.1,9.8,2.4,5.6,5.3),20)

myData <- data.frame (fReserve ,DivBoulders ,Roughness ,Biomass )

#glm
glm1 <- glm (Biomass ~ fReserve * Roughness + fReserve *  DivBoulders  ,
family= Gamma, data= myData)

#glmm:
library (nlme)

lme1 <- lme (Biomass  ~  fReserve , random= ~1 + fReserve | Roughness +
DivBoulders  , data=  myData ) # random intercept and slope
#Error in getGroups.data.frame(dataMix, groups) :
 # Invalid formula for groups
 # if I only use one random variable I have:
#Error in chol.default((value + t(value))/2) :
  #the leading minor of order 2 is not positive definite


lme2 <- lme( Biomass  ~  fReserve , random = ~1 | Roughness +
DivBoulders  ,data=myData) #random
intercept
#Error in getGroups.data.frame(dataMix, groups) :
 # Invalid formula for groups
# if I only use one random variable my result is fine!


lme3 <- lme (Biomass   ~  fReserve  , random=  ~ Roughness +   DivBoulders
| fReserve   , data= myData ) # from help (lme)
summary (lme3)
# I have a result. Is the model correct for what I want?
# BUT my fixed effect (reserve) does not have a p-value due to zero degrees
of freedom. However in glm1 it has.

Can you help me please?

Thanks a lot in advance,
Barbara

        [[alternative HTML version deleted]]