[R-br] Grafico de trajetórias - zap plot

Pedro Emmanuel Alvarenga Americano do Brasil emmanuel.brasil em gmail.com
Segunda Abril 1 23:28:28 BRT 2013


Amigos de R,

Estou fazendo umas análises iniciais de trajetórias de biomarcadores ao
longo de um tempo em que pacientes receberam um tratamento. Eu gostaria de
ter uma percepção do conjunto de trajetórias sem ajuste e procurei por ai
formas de faze-lo. Encontrei stats::interaction.plot e
epicalc::followup.plot. Mas por algum motivo desconhecido, essas duas
funções estavam truncando as trajetórias em 120 dias, mesmo solicitando
para mostrar as trajetórias em 180 dias que eu sei o que é o maximo de
tempo observado contido nos dados que tenho. Assim, depois de bater cabeça
um tanto investi um dia inteiro para fazer a função abaixo e compartilho
com voces.

Essa função é bem semelhante a followup.plot, mas, quando se deseja não
plotar todos os sujeitos para não poluir o gráfico, ao inves de selecionar
sujeitos aleatórios, eu coloquei um argumento chamado zap que seleciona n
sujeitos como valores iniciais, minimo, maximo e mediano.


# data = dataset in long format
# id = character indicating the name of an unique id variable for each
patient in data  (tested with numeric variable)
# time = a character indicating the name of a numeric variable representing
the time of the measure
# outcome = a character indicating the name of a numeric variable
representing the outcome
# plot.arg =  a list of arguments to pass to plot.default
# line.arg  = one or more  lists of arguments to pass to lines (one list to
each strata of by)
# zap = if auto then the n subjects with the lower, median, higher initial
outcome values will be ploted, anything different from "auto" will plot all
subjects
# n = number of subject to use in zap
# by = a character indicating the name of a factor/character variable
representing strata for ploting

# think about making conditions for zap such as only the median, only the
lowest or only the highest, the final highest, final lowest, final median
etc

zapplot <- function(data,id,time,outcome,
                    plot.arg=list(xlab='time',ylab=outcome),

line.arg=list(list(lty=2,col='blue'),list(lty=3,col='red')),
                    leg.args=list(x='top',inset=c(0,-.15),

horiz=TRUE,xjust=.5,col=c('blue','red'),lty=2:3,
                                  bty='n',xpd=T),
                    zap='auto',n=1,by=NULL){
#  if(!is.null(by)){
#    if(length(line.arg)!=nlevels(by)){
#      stop('If by is porvided, number of lists in line.arg and number of
levels in by should be the same! Try improving the lists in line.arg.')}
#  }
  data <- data[which(!is.na(data[,time]) & !is.na(data[,outcome])),]
  plot.arg <- c(list(x = data[,time],y = data[,outcome],type =
'n'),plot.arg)
  do.call(plot,plot.arg)
  if(!is.null(by)){
    if(!is.factor(data[,by])){
      data[,by] <- as.factor(data[,by])
    }
    tmp2 <- data
    for(j in 1:nlevels(data[,by])){ # i=2
      data <- tmp2
      data <- data[which(data[,by]==levels(data[,by])[j]),]
      tmp <- data[which(data[,time]==min(data[,time])),]
      tmp <- tmp[order(tmp[,outcome]),] # tmp[,c(id,time,outcome)]
      id.max <- tmp[(nrow(tmp)-(n-1)):nrow(tmp),id] #
      id.min <- tmp[(1:n),id] # data[which(data$id==id.min),outcome]
      tmp$med <- abs(tmp[,outcome]-median(tmp[,outcome],TRUE))
      tmp <- tmp[order(tmp$med),]
      id.med <- tmp[1:n,id]
      if(zap=='auto'){
        for(i in 1:length(id.min)){
          line.arg.min <- c(list(x =
data[which(data[,id]==id.min[i]),time],y =
data[which(data[,id]==id.min[i]),outcome]),line.arg[[j]])
          do.call(lines,line.arg.min)
        }
        for(i in 1:length(id.med)){
          line.arg.med <- c(list(x =
data[which(data[,id]==id.med[i]),time],y =
data[which(data[,id]==id.med[i]),outcome]),line.arg[[j]])
          do.call(lines,line.arg.med)
        }
        for(i in 1:length(id.max)){
          line.arg.max <- c(list(x =
data[which(data[,id]==id.max[i]),time],y =
data[which(data[,id]==id.max[i]),outcome]),line.arg[[j]])
          do.call(lines,line.arg.max)
        }
      }
      else{
        uni <- unique(data[,id])
        for(i in 1:length(uni)){
          line.arg.all <-
c(list(x=data[which(data[,id]==uni[i]),time],y=data[which(data[,id]==uni[i]),outcome]),line.arg[[j]])
          do.call(lines,line.arg.all)
        }
      }
    }
    # COLOCAR AQUI OS COMANDOS DE LEGENDA para identificar as diferentes
categorias do by
    leg.args <- c(legend=list(levels(data[,by])),leg.args)
    do.call(legend,leg.args)
  }
  else{
    tmp <- data[which(data[,time]==min(data[,time])),]
    tmp <- tmp[order(tmp[,outcome]),] # tmp[,c(id,time,outcome)]
    id.max <- tmp[(nrow(tmp)-(n-1)):nrow(tmp),id] #
    id.min <- tmp[(1:n),id] # data[which(data$id==id.min),outcome]
    tmp$med <- abs(tmp[,outcome]-median(tmp[,outcome],TRUE))
    tmp <- tmp[order(tmp$med),]
    id.med <- tmp[1:n,id]
    if(zap=='auto'){
      for(i in 1:length(id.min)){
        line.arg.min <- c(list(x = data[which(data[,id]==id.min[i]),time],y
= data[which(data[,id]==id.min[i]),outcome]),line.arg[[1]])
        do.call(lines,line.arg.min)
      }
      for(i in 1:length(id.med)){
        line.arg.med <- c(list(x = data[which(data[,id]==id.med[i]),time],y
= data[which(data[,id]==id.med[i]),outcome]),line.arg[[1]])
        do.call(lines,line.arg.med)
      }
      for(i in 1:length(id.max)){
        line.arg.max <- c(list(x = data[which(data[,id]==id.max[i]),time],y
= data[which(data[,id]==id.max[i]),outcome]),line.arg[[1]])
        do.call(lines,line.arg.max)
      }
    }
    else{
      uni <- unique(data[,id])
      for(i in 1:length(uni)){
        line.arg.all <-
c(list(x=data[which(data[,id]==uni[i]),time],y=data[which(data[,id]==uni[i]),outcome]),line.arg[[1]])
        do.call(lines,line.arg.all)
      }
    }
  }
}


# zapplot(b,'pront','dia','HLA',zap='all',n=2,by='CD4cod') #

Abraço forte e que a força esteja com voces,

Dr. Pedro Emmanuel A. A. do Brasil
Curriculum Lattes:  http://lattes.cnpq.br/6597654894290806
Instituto de Pesquisa Clínica Evandro Chagas
Fundação Oswaldo Cruz
Rio de Janeiro - Brasil
Av. Brasil 4365,
CEP 21040-360,
Tel 55 21 3865-9648
email: pedro.brasil em ipec.fiocruz.br
email: emmanuel.brasil em gmail.com

---Apoio aos softwares livres
www.zotero.org - gerenciamento de referências bibliográficas.
www.broffice.org ou www.libreoffice.org - textos, planilhas ou
apresentações.
www.epidata.dk - entrada de dados.
www.r-project.org - análise de dados.
www.ubuntu.com - sistema operacional
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20130401/0231f99b/attachment.html>


Mais detalhes sobre a lista de discussão R-br