0

我有来自几个人（“ ID” ）的单个变量（“ Count ”）与时间（“ Days ”）的纵向数据，以及为每个人生成带有单独的Count图的代码块。在某些情况下，一些数据点在变量“ Treat ”下标记为“1” ，并绘制为红色。其他的处理值 = 0。

我的问题：我需要 a) 为每个图添加一条回归线，仅使用Treat = 0 的数据点，并且 b) 需要在每个子图上打印线的斜率。我不知道该怎么做。到目前为止我的代码：

#plot longitudinal CD4 counts from multiple individuals in separate plots

#Clear work area
rm(list = ls(all = TRUE))

#Enter example data
tC <- textConnection("
ID  VisitDate   Count   Treat   Treatstarted
C0098   12-Feb-10   457 0   NA
C0098   2-Jul-10    467 0   NA
C0098   7-Oct-10    420 0   NA
C0098   3-Feb-11    357 0   NA
C0098   8-Jun-11    209 0   NA
C0098   9-Jul-11    223 0   NA
C0098   12-Oct-11   309 0   NA
C0110   23-Jun-10   629 0   30-Oct-10
C0110   30-Sep-10   461 0   30-Oct-10
C0110   15-Feb-11   270 1   30-Oct-10
C0110   22-Jun-11   236 1   30-Oct-10
C0151   2-Feb-10    199 0   NA
C0151   24-Mar-10   936 0   3-Apr-10
C0151   7-Jul-10    1147    1   3-Apr-10
C0151   9-Mar-11    1192    1   3-Apr-10
")
data1 <- read.table(header=TRUE, tC)
close.connection(tC)

# calculate elapsed time from first date, by ID
data1$VisitDate <- with(data1,as.Date(VisitDate,format="%d-%b-%y"))
data1$Days <- unlist(with(data1,tapply(VisitDate,ID,function(x){x-x[1]})))

#Define plot function
plot_one <- function(d){
 with(d, plot(Days, Count, t="n", tck=1, main=unique(d$ID), cex.main = 0.8, ylab = "", yaxt = 'n', xlab = "", xaxt="n",  xlim=c(0,1000), ylim=c(0,1200))) # set limits
    grid(lwd = 0.3, lty = 7)
    with(d[d$Treat == 0,], points(Days, Count, col = 1)) 
    with(d[d$Treat == 1,], points(Days, Count, col = 2))
}

#Create multiple plot figure
par(mfrow=c(8,8), oma = c(0.5,0.5,0.5,0.5), mar = c(0.5,0.5,0.5,0.5))
plyr::d_ply(data1, "ID", plot_one)

非常感谢的想法

4

1 回答 1

0

只需对您的功能稍作修改就可以plot_one解决问题。

plot_one <- function(d){
  with(d, plot(Days, Count, t="n", tck=1, main=unique(d$ID), cex.main = 0.8, ylab = "", yaxt = 'n', xlab = "", xaxt="n",  xlim=c(0,1000), ylim=c(0,1200))) # set limits
  grid(lwd = 0.3, lty = 7)
  with(d[d$Treat == 0,], points(Days, Count, col = 1)) 
  with(d[d$Treat == 1,], points(Days, Count, col = 2))
  mod = lm(Count ~ Days, data = d[d$Treat == 0,])
  abline(reg = mod)
  text(x=500, y=800, mod$coefficients[2])
}

这给出了以下图表（我更改了标准以使其更大一点）：

在此处输入图像描述

于 2013-11-15T15:31:49.810 回答