[R] Overlaying lattice graphs (continued)

[Deepayan Sarkar]

On 6/22/07, Sébastien <pomchip at free.fr> wrote:
> Hi Deepayan,
>
> The following code creates a dummy dataset which has the same similar as
> my usual datasets. I did not try to implement the changes proposed by
> Hadley, hoping that a solution can be found using the original dataset.
>
> ######### My code
>
> # Creating dataset
>
> nPts<-10            # number of time points
> nInd<-6              # number of individuals
> nModel<-3         # number of models
>
> TimePts<-rep(1:nPts,nInd*nModel)                                    #
> creates the "Time" column
> Coef<-rep(rnorm(6,0.1,0.01),each=nPts,nModel)             # Creates a
> vector of coefficients for generating the observations
> Obs<-10*exp(-Coef*TimePts)                                         #
> creates the observations
>
> for (i in 1:60){
> Pred[i]<-jitter(10*exp(-Coef[i]*TimePts[i]))
> Pred[i+60]<-jitter(5)
> Pred[i+120]<-jitter(10-Coef[i+120]*TimePts[i])
> }
>                   # creates the predicted values
>
> colPlot<-rep(1,nPts*nInd*nModel)
>     # creates the "Plot" column
> colModel<-gl(nModel,nPts*nInd,labels=c("A","B","C"))             #
> creates the "Model" column
> colID<-gl(nInd,nPts,nPts*nInd*nModel)
>       # creates the "ID" column
>
> mydata<-data.frame(colPlot,colModel,colID,TimePts,Obs,Pred)
>               # creates the dataset
> names(mydata)<-c("Plot","Model","Individuals","Time","Observed","Predicted")

The way you have structured your data makes no sense to me. In
particular, your 'Observed' data is the same set of 60 numbers
repeated 3 times, and this is not reflected in the data structure at
all. What would you want to happen if the numbers were not repeated?
Would you always plot the first 60, or would plot all of them?

If I understand what you are trying to do, this might be a more
transparent approach:


nPts<-10   # number of time points
nInd<-6    # number of individuals

TimePts <- rep(1:nPts, nInd)
Coef <- rep(rnorm(6,0.1,0.01), each = nPts)
Obs <- 10 * exp(-Coef * TimePts)
colID <- gl(nInd, nPts)

mydata <- data.frame(Time = TimePts, Observed = Obs, Individuals = colID)

fmA <- lm(Observed ~ Time, mydata)
fmB <- lm(Observed ~ poly(Time, 2), mydata)
fmC <- lm(Observed ~ poly(Time, 2) * Individuals, mydata)

mydata$PredA <- predict(fmA)
mydata$PredB <- predict(fmB)
mydata$PredC <- predict(fmC)

xyplot(Observed + PredA + PredB + PredC ~ Time | Individuals,
       data = mydata,
       type = c("p", "l", "l", "l"),
       distribute.type = TRUE)

-Deepayan