时间戳

首页 > 解决方案 > 如何对循环数据进行多路方差分析？

r - 如何对循环数据进行多路方差分析？

问题描述

我正在使用 R 包分析应变仪数据circular。我的数据是度数转换为弧度，我有四个自变量。其中一个 IV 有两个级别，但其他三个有 3-5 个级别。我已经很好地计算了数据的汇总统计数据（使用CircStats），但希望包括一个推理统计测试，如多路方差分析。该软件包circular具有这些数据的 ANOVA，aov.circular()但一次只针对一个 IV。我找不到使用循环数据的多路方差分析的任何等效项。有谁知道是否存在这样的测试，或者在任何现有的 R 包中是否可以进行这样的测试？

这是一个假设的数据框和汇总统计信息。

dfa <- data.frame(individual=c("C", "A", "B", "D", "D", "A", "A", "A"),
                  side=c("L", "L", "L", "R", "L", "L", "L", "R"),
                  type=c("ch", "ib", "ch", "ch", "cb", "ib", "pb", "pb"),
                  location=c("rear", "rear", "rear", "center", "center", "front", "front", "center"),
                  radian_degree=c(3.746309603, 3.875998892, 3.933472332, 3.253853158, 4.146252093, 4.522643462, 4.361775091, 3.759583424))

dfa %>% 
  group_by(individual, side, type, location) %>%
  summarise(circ.summary(radian_degree), est.kappa(radian_degree)) %>% 
  drop_na()

在您的带领下，我有以下...

dvm <- circular:::DvonmisesRad
stopifnot(all.equal(c(dvm(1,1,0.5,log=TRUE), dvm(2,3,0.5,log=TRUE)),
                    dvm(c(1,2), c(1,3), 0.5, log=TRUE)))

m_null <- mle2(radian_degree ~ dvm(exp(log_mu), exp(log_kappa)),
           parameters = list(log_mu ~ 1),
           data = df,
           start = list(log_mu = 0, log_kappa = 0))

m3_ILST <- update(m_null, parameters = list(log_mu ~ individual+location+side+type))
m3_ILS  <- update(m_null, parameters = list(log_mu ~ individual+location+side))
m3_IL  <- update(m_null, parameters = list(log_mu ~ individual+location))
m3_I  <- update(m_null, parameters = list(log_mu ~ individual))

model <- anova(m3_ILST, m3_ILS, m3_IL, m3_I, m_null)

但是，偏差都是相同的（即使包括了我的所有 9000 行数据）。此外，对于所有第 2-5 行的输出，Chisq 为 0，对于第 2-5 行的输出，Pr(>Chisq) 为 1。

标签： ranovamultivariate-testing

解决方案

原则上你可以这样做：

加载包：

library(tidyverse)
library(circular)
library(bbmle)

该circular::dvonmises()函数负责一系列参数检查和转换，但不允许我们计算向量化的对数似然mu（位置参数），因此我们将使用较低级别的内部版本：

dvm <- circular:::DvonmisesRad
## check vectorization
stopifnot(all.equal(c(dvm(1,1,0.5,log=TRUE), dvm(2,3,0.5,log=TRUE)),
                    dvm(c(1,2), c(1,3), 0.5, log=TRUE)))

拟合一个空模型。由于mu(location) 和kappa(scale) 参数都必须为非负数，因此我们将为两者使用日志链接。

m_null <- mle2(radian_degree ~ dvm(exp(log_mu), exp(log_kappa)),
           parameters = list(log_mu ~ 1),
           data = dfa,
           start = list(log_mu = 0, log_kappa = 0))

我们可以用不同的因素组合来拟合其他模型（+对于加法项*或:交互作用）

m3_ISL <- update(m1, parameters = list(log_mu ~ individual+side+location))
m3_IS  <- update(m1, parameters = list(log_mu ~ individual+side))
m3_I  <- update(m1, parameters = list(log_mu ~ individual))

然后我们可以anova()用来做序列似然比检验。

anova(m3_ISL, m3_IS, m3_I, m_null)
Likelihood Ratio Tests
Model 1: m3_ISL, radian_degree~dvm(exp(log_mu),exp(log_kappa)):
          log_mu~individual+side+location
Model 2: m3_IS, radian_degree~dvm(exp(log_mu),exp(log_kappa)):
          log_mu~individual+side
Model 3: m3_I, radian_degree~dvm(exp(log_mu),exp(log_kappa)): log_mu~individual
Model 4: m_null, radian_degree~dvm(exp(log_mu),exp(log_kappa)): log_mu~1
  Tot Df Deviance Chisq Df Pr(>Chisq)
1      8   29.406                    
2      6   29.406 0e+00  2     1.0000
3      5   29.406 1e-04  1     0.9939
4      2   29.406 1e-04  3     1.0000

这绝对不如罐装 ANOVA 表方便，但也不算太糟糕（我认为）。

这里的结果看起来很不稳定（偏差几乎相同），但我认为/希望这是因为数据集很小。这是一个表明它可以工作的示例：

set.seed(101)
dsim <- data.frame(f=factor(rep(1:2, each=20)),
                   radian_degree=unlist(lapply(c(pi/2, pi/4),
                                   rvonmises,
                                   n = 20,
                                   kappa = 20)))
ms_null <- update(m_null, data = dsim)
ms_f <- update(ms_null, parameters = list(log_mu ~ f))
anova(ms_null, ms_f)
Model 1: ms_null, radian_degree~dvm(exp(log_mu),exp(log_kappa)): log_mu~1
Model 2: ms_f, radian_degree~dvm(exp(log_mu),exp(log_kappa)): log_mu~f
  Tot Df Deviance  Chisq Df Pr(>Chisq)    
1      2   53.655                         
2      3    4.520 49.135  1  2.389e-12 ***

（我知道这个例子只是一个单向比较——你可以试着对照一下aov.circular()——但它表明这个建议并没有完全被误导......）

在多路设计中测试效果的一个合理但适度复杂（抱歉）的策略是：（1）通过测试有和没有它的模型来测试顶级（在您的情况下为 4 路）交互，即

parameters = list(log_mu ~ individual*side*location*type)
## vs.
parameters = list(log_mu ~ (individual+side+location+type)^3)

(2) 通过从模型中删除所有三向交互来测试每一个三向交互，例如

parameters = list(log_mu ~ (individual+side+location+type)^3)
## vs.
parameters = list(log_mu ~ (individual+side+location+type)^3-
                            individual:side:location)

测试 2-way 交互（仅包括不包含这些术语的更高级别的术语），例如

parameters = list(log_mu ~ (individual+side)^2 + location*type)
## vs.
parameters = list(log_mu ~ (individual+side)^2 + location + type)

等等。这种一般策略在?car::Anova

“type-II”和“type-III”的名称是从 SAS 借用的，但这里使用的定义与 SAS 使用的定义并不完全对应。Type-II测试是根据边缘性原则计算的，先测试每个term，除了忽略term的高阶亲属；所谓的 III 型测试违反了边际性，将模型中的每个项都测试在所有其他项之后。II 类检验的定义对应于 SAS 为方差分析模型生成的检验，其中所有预测变量都是因子，但不是更普遍的（即，当有定量预测变量时）。

推荐阅读