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@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