R epi friends,

Im conducting some initial longitudinal analysis of some biomarkers. I
would  like to have some idea of how the crude trajectories behave before
the attempt to run marginal models. I looked around and found the
stats::interaction.plot and epicalc::followup.plot. However, due to a
unknown reason, these functions were truncating the trajectories before the
maximum time. After a while banging my head around I decided to make my own
function, which I share with you below.

It is very similar to followup.plot, but in an attempt to make a
less polluted plot, instead of selecting randomly fewer subject it selects
n subjects closest to the initial median value, initial lower value and
initial maximum value.

In hope it may be useful to others too...

Kind regards,


# 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
# leg.args is a list of arguments to pass to legend, it will be used only
if by is different from NULL

# 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') #



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@ipec.fiocruz.br
email: emmanuel.brasil@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

	[[alternative HTML version deleted]]

