来自混合效应模型的 ggplot2 中的 ACF 图

ACF plot in ggplot2 from a mixed effects model

我正在尝试在 ggplot2 中执行自相关图。

library(nlme)

fm2 <- lme(distance ~ age + Sex, data = Orthodont, random = ~ 1)
plot(ACF(fm2,resType="normalized"),alpha=0.05)

通过以上函数得到的结果:

########################## IC ###########################

ic_alpha= function(alpha, acf_res){
  return(qnorm((1 + (1 - alpha))/2)/sqrt(acf_res$n.used))
}

#################### graphics ###########################
library(ggplot2)

ggplot_acf_pacf= function(res_, lag, label, alpha= 0.05){
  df_= with(res_, data.frame(lag, ACF))
  
  lim1= ic_alpha(alpha, res_)
  lim0= -lim1
  
  
  ggplot(data = df_, mapping = aes(x = lag, y = ACF)) +
    geom_hline(aes(yintercept = 0)) +
    geom_segment(mapping = aes(xend = lag, yend = 0)) +
    labs(y= label) +
    geom_hline(aes(yintercept = lim1), linetype = 2, color = 'blue') +
    geom_hline(aes(yintercept = lim0), linetype = 2, color = 'blue')
}
######################## result ########################
acf_ts = ggplot_acf_pacf(res_= ACF(fm2,resType="normalized"),
                         20,
                         label= "ACF")

但是,我遇到了以下错误:

Error in sqrt(acf_res$n.used) : 
  non-numeric argument to mathematical function

我打算得到的是这样的:

ACF 生成的对象 没有 名为 n.used 的成员。它有一个名为 n.used 属性 。所以你的 ic_alpha 函数应该是:

ic_alpha <- function(alpha, acf_res) {
  return(qnorm((1 + (1 - alpha)) / 2) / sqrt(attr(acf_res, "n.used")))
}

另一个问题是,由于 ic_alpha returns 一个向量,你不会有一对显着性线,而是每个滞后一对,这看起来很乱。相反,模拟基本的 R 绘图方法,我们可以使用 geom_line 来获得单个曲线对。

ggplot_acf_pacf <- function(res_, lag, label, alpha = 0.05) {
  
  df_ <- with(res_, data.frame(lag, ACF))
  
  lim1 <- ic_alpha(alpha, res_)
  lim0 <- -lim1
  
  ggplot(data = df_, mapping = aes(x = lag, y = ACF)) +
    geom_hline(aes(yintercept = 0)) +
    geom_segment(mapping = aes(xend = lag, yend = 0)) +
    labs(y= label) +
    geom_line(aes(y = lim1), linetype = 2, color = 'blue') +
    geom_line(aes(y = lim0), linetype = 2, color = 'blue') +
    theme_gray(base_size = 16)
}

这导致:

ggplot_acf_pacf(res_ = ACF(fm2, resType = "normalized"), 20, label = "ACF")