各位网友：
我在读一篇论文（附件1），但是在如何产生文中的Figure 1是遇到了难题。我按照文中第2部分的产生数据的方法产生了所需数据，但是利用附件2 给出的R代码并不能产生图1， 所以请求各位指点迷津。
产生数据代码：
##set the seed to make the result is reproducible;
set.seed(100);
##specify the number of the subjects in the trial;
n<-100
##specify the number of visits planned;
v<-5
##sample the measures as the complete data;
x<-matrix(rlnorm(n*v,0,1), nrow=n);#sample the complete data; n subjects v visits planned

##rename the complete data so that it can be handled later;
y<-x;
##specify the drop out timepoints for each subject independently (sample should with replacement);
drop<-sample(1:v, n, replace=T)
#rename the drop timepoints;
d<-drop;
#change the data to make the data meet the drop meaning;
for(subj in 1:n){
while(d[subj]<v){
y[subj, d[subj]+1]=NA
d[subj]=d[subj]+1 }
}
##disply the final result;
y

#calculate the mean measurement by drop out time and meansureable timepoint;
m<-matrix(NA, nrow=v, ncol=v);
#redefine the matrix of the data Xit*I(Di=d)
ind<-matrix(, nrow=1,ncol=v)
#Get the final data as Table 1;

for (dr in 1:v){
for(tr in 1:v){
if (length(which(drop==dr))>0){ind<-which(drop==dr)
m[v-dr+1,tr]<-mean(y[ind,tr])}
}
}
##Get the final data to be used Plot;
m
附件中的产生图1的代码：
require(plotrix)
first.na=function(x){
if (sum(is.na(x))==0) return(length(x)+1)
else return(min(which(is.na(x))))
}
last.non.na=function(x){
return(which.max(!is.na(x)==T))
}
triPlot=function(dsn,incl.leg=F,legloc=rep(0,4),cex=1,maintitle=""){
tt=ncol(dsn)
mm=apply(dsn,1,last.non.na)
timemeans=matrix(NA,nrow=tt,ncol=tt)
for(j in 1:tt){
timemeans[j,]=apply(dsn[which(mm==j),],2,mean,na.rm=T)
}
grcolor=grey(90:1 / 100)

par(mar=c(4,4,6,2)+0.1,mgp=c(3,0.5,0))
image(t(timemeans),x=1:tt,y=1:tt,col=grcolor,xaxt="n",yaxt="n",bty="n",
xlab="",ylab="Dropout Time",main=maintitle)
abline(v=(1:tt)+0.5,col="white")
abline(h=(1:tt)+0.5,col="white")
axis(2,lwd=0)
axis(3,lwd=0)
mtext("Measurement Time",side=1,line=1)
if(incl.leg)
color.legend(legloc[1],legloc[2],legloc[3],legloc[4],
legend=pretty(timemeans,n=2,min.n=1),
rect.col=grcolor,gradient="x",cex=cex)
}

运行 triPlot(m) 只能产生最后的 图， 请指导。谢谢。