時系列データの描画

同僚から、時系列データの解析を頼まれた。線形混合モデルを勧めたが、とりあえず図示することが大事である。

どんなに良いモデル選択法を用いても、やはり人間の直感/見た目は重要である。

しかーし、良い関数を知らないだけかもしれないが欠損値が含まれる時系列データをうまく図示することができなかった。

仕方ないので、Rの勉強を兼ねて自作関数を作った。しかし、はじめは意図とは違う感じになり、相当に苦労した。目的は達成したものの我ながら、冗長な感じで今ひとつだ。先輩方でもっとよい実装法があればご教授いただければ嬉しいです。

欠損値がなければ簡単だったが、まだらに欠損値があっても図示出来るようにした。欠損値があれば、そこは飛ばして、線でつなぐようにする。 long型用とwide型用を作る。（reshapeすれば良いんだけどね。）

まずはwide型の時系列データを図示する。

set.seed(1)
n <- 40
a1 <- a2 <- a3 <- a4 <- numeric(n)
id <- time <- y <- r <- numeric(n * 4)
a1 <- rnorm(n, 0, 10)
a2 <- rnorm(n, -5, 10)
a3 <- rnorm(n, 10, 10)
a4 <- rnorm(n, 5, 10)
id <- rep(1:n, 4)
time <- as.numeric(gl(4, n))
y <- c(a1, a2, a3, a4)
long1 <- data.frame(id, time, y)
r <- runif(n * 4)

for (i in 1:(4 * n)) {
    if (r[i] >= 0.66) {
        long1[i, ] <- NA
    }
}
long1 <- subset(long1, complete.cases(long1))

wide1 <- reshape(long1, timevar = "time", idvar = "id", direction = "wide")



plot.by.time <- function(var, data) {
    # varに時系列データにしたい列を順番にすべて入れる。dataはデータフレーム
    boxplot(data[, var])
    cols = rainbow(nrow(data), start = 0.5, end = 0.1)  #色の設定を行う


    for (l in 1:nrow(data)) {
        for (m in 1:length(var)) {
            points(m, data[l, var[m]], col = cols[l], pch = 16)
        }
    }



    for (j in seq(nrow(data))) {
        for (i in (seq(length(var) - 1))) {
            comp <- 0
            if (is.na(data[j, var[i]])) {
                next
            } else {
                for (k in (2:length(var))) {
                  if (i >= k) {
                    next
                  } else if (is.na(data[j, var[k]])) {
                    next
                  } else if (comp == 1) {
                    break
                  } else {
                    segments(i, data[j, var[i]], k, data[j, var[k]], lty = 2, 
                      col = cols[j])
                    comp <- 1
                  }
                }
            }
        }
    }
}


plot.by.time(c("y.1", "y.2", "y.3", "y.4"), wide1)


## long形式####


plot.by.time.long <- function(out, time, time_0 = 0, id, data) {
    # outcome, data, id, timeを表す変数 time0は原点,timeは必ず連続数値

    data <- data[order(data[, id], data[, time]), ]  #id,time順に並び替える
    ud <- unique(data[, id])  #id の種類の数

    cols = rainbow(length(unique(data[, id])), start = 0.5, end = 0.1)  #色の設定を行う
    # boxplot(data[,out]~data[,time])
    plot(2, 0, xlim = c(min(data[, time] - time_0, na.rm = TRUE), max(data[, 
        time] - time_0, na.rm = TRUE)), ylim = c(min(data[, out], na.rm = TRUE), 
        max(data[, out], na.rm = TRUE)), type = "n")
    for (l in 1:nrow(data)) {
        points(data[l, time] - time_0, data[l, out], col = cols[which(ud == 
            data[l, id])], pch = 16, xlim = c(min(data[, time] - time_0, na.rm = TRUE), 
            max(data[, time] - time_0, na.rm = TRUE)), ylim = c(min(data[, out], 
            na.rm = TRUE), max(data[, out], na.rm = TRUE)))
    }

    for (j in 1:(nrow(data) - 1)) {
        if (data[j, id] == data[j + 1, id]) {
            segments(data[j, time] - time_0, data[j, out], data[j + 1, time] - 
                time_0, data[j + 1, out], lty = 2, cols[which(ud == data[j, 
                id])], xlim = c(min(data[, time] - time_0, na.rm = TRUE), max(data[, 
                time] - time_0, na.rm = TRUE)), ylim = c(min(data[, out], na.rm = TRUE), 
                max(data[, out], na.rm = TRUE)))
        }
    }

}


plot.by.time.long(out = "y", time = "time", time_0 = 0, id = "id", data = long1)

plot of chunk unnamed-chunk-1

なんとかできた

r-statistics-fanの日記

統計好き人間の覚書のようなもの

時系列データの描画

時系列データの描画