在ggplot2中，使用"lomb"包绘制的图表

Question

我正在使用“lomb”包来计算Lomb-Scargle周期图，这是一种分析生物时间序列数据的方法。如果您要求它这样做，该包确实会创建一个图形。然而，与ggplot2图形相比，这些图形不太好看。因此，我想用ggplot绘制结果。但是，我不知道如何访问绘制曲线的函数...

这是一个绘图的示例代码:

TempDiff <- runif(4033, 3.0, 18) % just generate random numbers
Time2 <- seq(1,4033) % Time vector
Rand.LombScargle <- randlsp(repeats=10, TempDiff, times = Time2, from = 12, to = 36, 
                        type = c("period"), ofac = 10, alpha = 0.01, plot = T, 
                        trace = T, xlab="period", main = "Lomb-Scargle Periodogram")

我也尝试着查找函数randlsp本身的信息，但是在那里我并没有找到任何有用的东西...

 getAnywhere(randlsp)
 A single object matching ‘randlsp’ was found
 It was found in the following places
   package:lomb
   namespace:lomb
 with value

 function (repeats = 1000, x, times = NULL, from = NULL, to = NULL, 
 type = c("frequency", "period"), ofac = 1, alpha = 0.01, 
 plot = TRUE, trace = TRUE, ...) 
  {
    if (is.ts(x)) {
    x = as.vector(x)
  }
 if (!is.vector(x)) {
    times <- x[, 1]
    x <- x[, 2]
}
 if (plot == TRUE) {
    op <- par(mfrow = c(2, 1))
}
 realres <- lsp(x, times, from, to, type, ofac, alpha, plot = plot, 
    ...)
 realpeak <- realres$peak
 pks <- NULL
 if (trace == TRUE) 
    cat("Repeats: ")
 for (i in 1:repeats) {
    randx <- sample(x, length(x))
    randres <- lsp(randx, times, from, to, type, ofac, alpha, 
        plot = F)
    pks <- c(pks, randres$peak)
    if (trace == TRUE) {
        if (i/10 == floor(i/10)) 
            cat(i, " ")
    }
}
if (trace == TRUE) 
    cat("\n")
prop <- length(which(pks >= realpeak))
p.value <- prop/repeats
if (plot == TRUE) {
    mx = max(c(pks, realpeak)) * 1.25
    hist(pks, xlab = "Peak Amplitude", xlim = c(0, mx), main = paste("P-value: ", 
        p.value))
    abline(v = realpeak)
    par(op)
}
res = realres[-(8:9)]
res = res[-length(res)]
res$random.peaks = pks
res$repeats = repeats
res$p.value = p.value
class(res) = "randlsp"
return(invisible(res))

任何想法都将受到赞赏！

最好的祝福， 克里斯汀

附：带有真实数据的图表示例。

- Christine Blume

---

说实话，我认为Lomb-Scargle周期图最初是在天体物理学中开发的（请参阅Press等人的《Numerical Recipes》）。 - Ben Bolker

Answer 1

将任何返回对象中的 ggplot 图形提取出来的关键是将所需数据转换为某种形式的 data.frame。为此，您可以查看返回值的对象类型，并查看可以立即提取到 data.frame 中的数据类型。

str(Rand.LombScargle) # get the data type and structure of the returned value

    List of 12
 $ scanned     : num [1:2241] 12 12 12 12 12 ...
 $ power       : num [1:2241] 0.759 0.645 0.498 0.341 0.198 ...
 $ data        : chr [1:2] "times" "x"
 $ n           : int 4033
 $ type        : chr "period"
 $ ofac        : num 10
 $ n.out       : int 2241
 $ peak        : num 7.25
 $ peak.at     : num [1:2] 24.6908 0.0405
 $ random.peaks: num [1:10] 4.99 9.82 7.03 7.41 5.91 ...
 $ repeats     : num 10
 $ p.value     : num 0.3
 - attr(*, "class")= chr "randlsp"

在randlsp的情况下，它是一个列表，通常是从统计函数返回的。大部分信息也可以从?randlsp中获取。

看起来Rand.LombScargle$scanned和Rand.LombScargle$power包含了第一个图所需的大部分内容：

Periodogram上还有一条水平线，但它与randlsp返回的任何内容都不对应。查看您提供的源代码，似乎Periodogram实际上是由lsp()生成的。

LombScargle <- lsp( TempDiff, times = Time2, from = 12, to = 36, 
    type = c("period"), ofac = 10, alpha = 0.01, plot = F)

str(LombScargle)

List of 12
 $ scanned  : num [1:2241] 12 12 12 12 12 ...
 $ power    : num [1:2241] 0.759 0.645 0.498 0.341 0.198 ...
 $ data     : chr [1:2] "Time2" "TempDiff"
 $ n        : int 4033
 $ type     : chr "period"
 $ ofac     : num 10
 $ n.out    : int 2241
 $ alpha    : num 0.01
 $ sig.level: num 10.7
 $ peak     : num 7.25
 $ peak.at  : num [1:2] 24.6908 0.0405
 $ p.value  : num 0.274
 - attr(*, "class")= chr "lsp"

根据这些数据，我猜测这条线表示的是显著性水平LombScargle$sig.level

将这些内容整合起来，我们可以从lsp中创建数据，以传递给ggplot：

lomb.df <- data.frame(period=LombScargle$scanned, power=LombScargle$power)

# use the data frame to set up the line plot
g <- ggplot(lomb.df, aes(period, power)) + geom_line() + 
       labs(y="normalised power", title="Lomb-Scargle Periodogram")

# add the sig.level horizontal line
g + geom_hline(yintercept=LombScargle$sig.level, linetype="dashed")

对于直方图，它看起来是基于randlsp中的向量Rand.LombScargle$random.peaks：

rpeaks.df <- data.frame(peaks=Rand.LombScargle$random.peaks)

ggplot(rpeaks.df, aes(peaks)) + 
       geom_histogram(binwidth=1, fill="white", colour="black") + 
       geom_vline(xintercept=Rand.LombScargle$peak, linetype="dashed") +
       xlim(c(0,12)) +  
       labs(title=paste0("P-value: ", Rand.LombScargle$p.value), 
           x="Peak Amplitude", 
           y="Frequency")

尝试调整这些图表，使其符合您的口味。

- MattLBeck

---

非常完美，非常感谢！您的解决方案真正帮助我理解了如何解决这个问题，并为我提供了能力，使我能够自己处理这种问题。 - Christine Blume