How to produce the Figure

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

各位网友：
我在读一篇论文（附件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) 只能产生最后的 图， 请指导。谢谢。

supp app.pdf (31.61 KB)

2013-8-19 15:16:29 上传

Visualizing Longitudinal Data with Drop Outs.pdf (273.32 KB)

2013-8-19 15:16:16 上传

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：Produce produc Figure prod The produce

Nobody can help me do this out?

can you help me?

问题出在 last.non.na 这个函数上：

1. last.non.na=function(x){
   
2.   return(which.max(!is.na(x)==T))
   
3. }

复制代码

该函数是为了得到一个对象重复观测中最后一次（缺失前）的观测位置，应该和函数 first.na() 类似改为：

1. last.non.na=function(x){
   
2.   return(max(which(!is.na(x)==T)))
   
3. }

复制代码

然后就是 triPlot() 这个函数，里面的 timemeans 矩阵似乎就是前面数据生成最后得到的m，这样把 y 当成 triPlot 的参数就可以得到下面的图：


另外你用m做参数也应该得到一样的图，但是要把 triPlot() 这个函数中
timemeans[j,]=apply(dsn[which(mm==j),],2,mean,na.rm=T) 这句完善下：

1. timemeans[j,]=apply(dsn[which(mm==j), ,drop=FALSE],2,mean,na.rm=T)

复制代码

加上红色的,drop=FALSE是为了防止1Xn矩阵降维成向量。

这样运行 triPlot(y) 和 triPlot(m) 得到的图就是一样的啦

树袋熊2 发表于 2013-8-21 20:46
问题出在 last.non.na 这个函数上：

该函数是为了得到一个对象重复观测中最后一次（缺失前）的观测位置， ...

Thanks very much, really appreicated for your kind help