[R-sig-ME] Odd plot.lme() behavior

[David Daniel]

Hi all,

I've come across what appears to me to be strange behavior in the  
plot.lme() method.  In the example code below, the first plot() call  
gives a plot with no points plotted.  But the second plot() call gives  
a plot as expected (with points).  The only difference is that in the  
second case, lme() uses a data set not containing the variable, YES,  
but this variable is not referenced in any of the functions used in  
the code.

Even more odd is that if a seed of 1 is used in the call to  
set.seed(), the second call to plot() gives an error message not given  
when the seed is 2 (the seed only affects which values of variables  
YES and NO are NA's, and these variables are not used):

|> plot(lmeB,  Group ~ resid(.,  type="p"))
|  Error in `[<-.data.frame`(`*tmp*`, , ".y", value = c(1, 1, 1, 1, 1,  
1,  :
|  replacement has 36 rows, data has 80

But if *both* variables YES and NO are removed from the data set used,  
the plot works fine with a seed of 1.

Have others seen this?  Is this expected behavior?

I'm on a MacBook Pro 2.16 GHz Core 2 Duo, with OS 10.5.5, R version  
2.7.0, package nlme version 3.1-89, and package lattice version  
0.17-8.  The sessionInfo() output is below.

Thanks,
-David

#---------------------------------------
library(lattice);
library(nlme);

# Save the user's random seed
save.seed <- .Random.seed


#---------------------------------------
set.seed(2, kind = NULL)

n <- 10	# Number of subjects in each Group
k <- 4	# Number of phases

Group <- factor(c(rep("A",n*k), rep("B",n*k)))
m <- length(Group)/k		# total subjects

Gender <- factor( rep(c(rep("M", k), rep("F", k)), n) )
ID <- rep(1:m, rep(k, m))
phase <- rep(1:k, m)
YES <- rep(ceiling(12*runif(m)), rep(k, m))
NO <- YES
YES[YES > 6] <- NA
NO[NO <= 6] <- NA
PB <- rnorm(m*k) + phase + as.numeric(Group)
plotData <- data.frame(cbind(ID, PB, Group, Gender, phase, YES, NO))

#---------------------------------------
# Removing either YES or NO from the data.frame gives the correct plot
#---------------------------------------
plotData2 <- subset(plotData, select=c(ID, PB, Group, phase, YES, NO))
lmeA <- lme(PB ~ Group + phase,
					random= list(ID= pdDiag(form= ~phase)),
					data=plotData2,
					na.action= na.omit)
plot(lmeA,  Group ~ resid(.,  type="p"))

#---------------------------------------
# Removing either YES or NO from the data.frame gives the correct plot
#---------------------------------------
plotData3 <- subset(plotData, select=c(ID, PB, Group, Gender, phase,  
YES))
lmeB <- lme(PB ~ Group + phase,
					random= list(ID= pdDiag(form= ~phase)),
					data=plotData3,
					na.action= na.omit)
plot(lmeB,  Group ~ resid(.,  type="p"))

#---------------------------------------
# Restore the user's random seed
.Random.seed <- save.seed
#---------------------------------------



 > sessionInfo()
R version 2.7.0 (2008-04-22)
i386-apple-darwin8.10.1

locale:
en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

loaded via a namespace (and not attached):
[1] Matrix_0.999375-9 grid_2.7.0        lattice_0.17-8     
lme4_0.999375-20
[5] nlme_3.1-89

----------------------------------
David Daniel
Associate Professor
University Statistics Center
New Mexico State University

ddaniel at nmsu.edu