在 2-y 轴散点图中添加单侧重复测量检验的中值趋势线和 p 值 [R]

Question

加载示例数据框

df <- structure(list(ID = c(1,1,1,2,2,2,3,3,3), 
time = c(0L,1L,2L,0L,1L,2L,0L,1L,2L),
M1a = c(0, 0.2, 0.3, 0, 1.5, 2.9,0, 2.4, 3.9), 
M2a = c(0, 0.4, 0.6,0,0.9, 0.9,0,0.5, 0.7), 
M3a = c(0,0.3, 0.4, 0, 0.6, 0.9,0, 0.5, 0.8), 
M4a = c(0,0.6, 0.6,0, 0.4, 0.6,0, 0.2, 0.9), 
M1b = c(0L, 200L, 300L,0L, 300L, 900L,0L, 900L, 1000L), 
M2b = c(0L, 400L, 600L,0L, 600L, 900L,0L, 600L, 1000L), 
M3b = c(0L, 300L, 400L,0L, 200L, 800L,0L, 200L, 900L), 
M4b = c(0L, 600L, 600L,0L, 800L, 1000L,0L, 400L, 1100L)), 
.Names = c("ID", "time", "M1a", "M2a", "M3a", "M4a","M1b", "M2b","M3b", "M4b"), class = "data.frame", row.names = c(NA, -9L))

现在绘制两个 y 轴散点图

par(mar=c(5,4,4,5)+.1)
plot(df$time,df$M1a,type="p",col="red", main="M1", cex=0.5, cex.main=2, cex.lab=1.0, cex.axis=0.7)
par(new = TRUE)
plot(df$time,df$M1b,type="p",col="blue",xaxt="n",yaxt="n",xlab="",ylab="")
mtext("Relative change (%)",side=4,line=3)
axis(4)
legend("topleft",col=c("red","blue"),lty=1,legend=c("Absolute Change","Relative Change"))

我被什么困住了？

1.中值趋势线

我能够添加回归线，但我想要一条连接三个时间点的 M1a 和 M1b 中位数的中位数趋势线。

2.在图中添加p值，重复单向方差分析

fit1=aov(df$M1a~df$time + Error(ID/time),na.action=na.exclude,data=df);
sig1= summary(fit1)$"Error: Within"$"Pr(>F)"
if (sig<0.001) star='**' else if (sig>=0.001&sig<0.05) star='*' else star='';
if (sig1<0.001) star='**' else star='';

我计划在我的 2-y 轴图中添加使用上述代码来添加 p 值。在这里，我将 sig1 设为 NULL，但是，sig1 应该打印出 0.153。

如果结果显着，最终结果应在情节 (M1) 的主标题上包含 * 标记。

有什么建议吗？提前致谢！

Answer 1

首先要回答 #2，需要查看 summary.aov 对象的内部结构：

dput(summary(fit1))
structure(list(`Error: ID` = structure(list(structure(list(Df = 1, 
    `Sum Sq` = 5.60666666666667, `Mean Sq` = 5.60666666666667, 
    `F value` = NA_real_, `Pr(>F)` = NA_real_), .Names = c("Df", 
"Sum Sq", "Mean Sq", "F value", "Pr(>F)"), class = c("anova", 
"data.frame"), row.names = "Residuals")), class = c("summary.aov", 
"listof")), `Error: ID:time` = structure(list(structure(list(
    Df = 1, `Sum Sq` = 11.3157142857143, `Mean Sq` = 11.3157142857143), .Names = c("Df", 
"Sum Sq", "Mean Sq"), class = c("anova", "data.frame"), row.names = "df$time")), class = c("summary.aov", 
"listof")), `Error: Within` = structure(list(structure(list(Df = c(1, 
5), `Sum Sq` = c(0.325952380952381, 0.573888888888889), `Mean Sq` = c(0.325952380952381, 
0.114777777777778), `F value` = c(2.83985617480293, NA), `Pr(>F)` = c(0.152766396924706, 
NA)), .Names = c("Df", "Sum Sq", "Mean Sq", "F value", "Pr(>F)"
), class = c("anova", "data.frame"), row.names = c("df$time  ", 
"Residuals"))), class = c("summary.aov", "listof"))), .Names = c("Error: ID", 
"Error: ID:time", "Error: Within"), class = "summary.aovlist")

请注意，summary(fit1)$"Error: Within" 中的值实际上被埋得更深一层（并且没有名称，因此需要数字索引。这样做：

summary(fit1)$"Error: Within"[[1]]$`Pr(>F)`
[1] 0.1527664        NA

现在看看我是否能弄清楚二-0 坐标的情节问题。很确定需要在par(new=TRUE) 操作之前进行任何中值绘图，因为这会根据新数据更改用户坐标系。

为您的情节添加具有提取价值的标题，并由@VincentBonhomme 的有用评论增强：

 plot(df$time,df$M1a,type="p",col="red", cex=0.5, cex.main=2, cex.lab=1.0, cex.axis=0.7)
lines(unique(df$time), 
        tapply(df$M1a, df$time, median)) 
par(new = TRUE) 
plot( df$time, df$M1b,type="p", col="blue", xaxt="n", yaxt="n", xlab="",ylab="") 
lines(unique(df$time), 
      tapply(df$M1b, df$time, median))
mtext("Relative change (%)",side=4,line=3)
axis(4)
legend("topleft",col=c("red","blue"), lty=1,legend=c("Absolute Change","Relative Change"))

title(main=bquote("P-value for M1 (absolute scale)"== 
                       .(round(summary(fit1)$"Error: Within"[[1]]$`Pr(>F)`, 3) ) )  )