R help - Sep 2003 - lme in R and Splus

Good Day,

Included below is some code to generate data and to fit a mixed effects
model to this fake data.  The code works as expected when I call the
function "lme" in Splus but not in R.  

The error message from calling lme in R is: 
"Error in getGroups.data.frame(dataMix, groups) :
        Invalid formula for groups"

I installed the nlme package for R around 20 August 2003. 

Thanks in advance.

System information:

Splus:
Version 6.1.2 Release 2 for Sun SPARC, SunOS 5.6 : 2002

R:
platform i686-pc-linux-gnu
arch     i686
os       linux-gnu
system   i686, linux-gnu
status
major    1
minor    7.1
year     2003
month    06
day      16
language R



############## BEGINNING OF CODE ###########################
# a fake dataset to make the bumps with
nn <- 30  # of data points
mm <- 7   # number of support sites for x(s)
# create sites s
ss <- seq(1,10,length=nn)
# create the data y
e1 <- rnorm(nn,sd=0.1)
e2 <- cos(ss/10*2*pi*4)*.2
yy <- sin(ss/10*2*pi)+e2+e1
plot(ss,yy)

# locations of support points
ww <- seq(1-2,10+2,length=mm)
# width of kernel
sdkern <- 2

# create the matrix KK
KK <- matrix(NA,ncol=mm,nrow=nn)
for(ii in 1:mm){
KK[,ii] <- dnorm(ss,mean=ww[ii],sd=sdkern)
}

# create a dataframe to hold the data
df1 <- data.frame(y=yy,K=KK,sub=1)
df1$sub <- as.factor(df1$sub)

# now fit a mixed model using lme
a1 <- lme(fixed= y ~ 1,
          random= pdIdent(~KK-1),
          data=df1,na.action=na.omit)

# obtain and plot the fitted values
a1p <- as.vector(predict(a1,df1))
lines(ss,a1p,lty=1)

##################### END OF CODE ######################################3

-- 
*********************************************************************
| Michael Fugate                         Temp Phone: (505) 665-1817 |
| Statistical Sciences Group, D-1                                   |
| Los Alamos National Laboratory         email: fugate at lanl.gov     |
| Los Alamos, NM 87545                                              |
| Mail Stop: F600                                                   |

Michael Fugate <fugate at lanl.gov> writes:
> ############## BEGINNING OF CODE ###########################
> # a fake dataset to make the bumps with
> nn <- 30  # of data points
> mm <- 7   # number of support sites for x(s)
> # create sites s
> ss <- seq(1,10,length=nn)
> # create the data y
> e1 <- rnorm(nn,sd=0.1)
> e2 <- cos(ss/10*2*pi*4)*.2
> yy <- sin(ss/10*2*pi)+e2+e1
> plot(ss,yy)
> 
> # locations of support points
> ww <- seq(1-2,10+2,length=mm)
> # width of kernel
> sdkern <- 2
> 
> # create the matrix KK
> KK <- matrix(NA,ncol=mm,nrow=nn)
> for(ii in 1:mm){
> KK[,ii] <- dnorm(ss,mean=ww[ii],sd=sdkern)
> }
> 
> # create a dataframe to hold the data
> df1 <- data.frame(y=yy,K=KK,sub=1)
> df1$sub <- as.factor(df1$sub)
> 
> # now fit a mixed model using lme
> a1 <- lme(fixed= y ~ 1,
>           random= pdIdent(~KK-1),
>           data=df1,na.action=na.omit)
You don't have a grouping factor in the random specification and I
can't tell from the simulation what you would expect the groups to be.
> # obtain and plot the fitted values
> a1p <- as.vector(predict(a1,df1))
> lines(ss,a1p,lty=1)
> 
> ##################### END OF CODE ######################################3
> 
> -- 
> *********************************************************************
> | Michael Fugate                         Temp Phone: (505) 665-1817 |
> | Statistical Sciences Group, D-1                                   |
> | Los Alamos National Laboratory         email: fugate at lanl.gov     |
> | Los Alamos, NM 87545                                              |
> | Mail Stop: F600                                                   |
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
-- 
Douglas Bates                            bates at stat.wisc.edu
Statistics Department                    608/262-2598
University of Wisconsin - Madison        http://www.stat.wisc.edu/~bates/

Michael Fugate <fugate at lanl.gov> writes:
> ############## BEGINNING OF CODE ###########################
> # a fake dataset to make the bumps with
> nn <- 30  # of data points
> mm <- 7   # number of support sites for x(s)
> # create sites s
> ss <- seq(1,10,length=nn)
> # create the data y
> e1 <- rnorm(nn,sd=0.1)
> e2 <- cos(ss/10*2*pi*4)*.2
> yy <- sin(ss/10*2*pi)+e2+e1
> plot(ss,yy)
> 
> # locations of support points
> ww <- seq(1-2,10+2,length=mm)
> # width of kernel
> sdkern <- 2
> 
> # create the matrix KK
> KK <- matrix(NA,ncol=mm,nrow=nn)
> for(ii in 1:mm){
> KK[,ii] <- dnorm(ss,mean=ww[ii],sd=sdkern)
> }
> 
> # create a dataframe to hold the data
> df1 <- data.frame(y=yy,K=KK,sub=1)
> df1$sub <- as.factor(df1$sub)
> 
> # now fit a mixed model using lme
> a1 <- lme(fixed= y ~ 1,
>           random= pdIdent(~KK-1),
>           data=df1,na.action=na.omit)
> 
> # obtain and plot the fitted values
> a1p <- as.vector(predict(a1,df1))
> lines(ss,a1p,lty=1)
lme in S-PLUS is older than the one in R, and some things changed. I
think you want

df1 <- data.frame(y=yy,K=I(KK),sub=1)
a1 <- lme(fixed= y ~ 1,
          random= list(sub=pdIdent(~K-1)),
          data=df1,na.action=na.omit)
lines(ss,predict(a1,df1,1))

(Apparently you can't do a level-0 prediction in a model with only an
intercept, which looks like a bit of a bug. Of course, that is just
the intercept for all observations, but...)

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907