标签: anova

我正在处理一个不平衡的设计/样本并且最初学到了aov().我现在知道,对于我的ANOVA测试,我需要使用III型平方和,这涉及使用拟合lm()而不是使用aov().

问题是使用事后测试(特别是Tukey的HSD)lm().我所做的所有研究都表示simint在multcomp包中使用会起作用,但现在已经更新了该命令似乎无法使用.它似乎也依赖于aov()计算.

基本上我为R找到的所有Tukey HSD测试都假定你aov()用于比较而不是lm().为了获得我必须使用的不平衡设计所需的III型平方和:

mod<-lm(Snavg~StudentEthnicity*StudentGender)

Anova(mod, type="III")

如何使用我的mod使用Tukey HSD测试lm()？或者相反,使用Type III计算我的ANOVA并仍然能够进行Tukey HSD测试？

谢谢!

试图学习R.来自旧统计数据的问题想要知道不同建筑工地的休息时间是否有差异.麻烦的是,文本决定每个站点雇用不同数量的工人.所以,我被困在使用不等样本大小的ANOVA寻求帮助.

site1 <- c(34,25,27,31,26,34,21)
site2 <- c(33,35,31,31,42,33)
site3 <- c(17,30,30,26,32,28,26,29)
site4 <- c(28,33,31,27,32,33,40)

我试图用SKLearn做一个LR,用于一个相当大的数据集,其中有600个虚拟数据集,只有很少的区间变量(我的数据集中有300 K行),结果混淆矩阵看起来很可疑.我想检查返回系数和ANOVA的重要性,但我找不到如何访问它.有可能吗？对于包含大量虚拟变量的数据,最佳策略是什么？非常感谢!

是aov适合不平衡数据集.根据帮助...provides a wrapper to lm for fitting linear models to balanced or unbalanced experimental designs.但后来就说了aov is designed for balanced designs, and the results can be hard to interpret without balance.

我应该如何在R中的非平衡数据集上执行双向anova？

我想重现SAS输出的I型和III型平方和的不同结果(使用时proc glm).我记得我们type III sum of squares用于不平衡的数据集.

先感谢您.

我目前正在测试是否应该在我的lmer模型中包含某些随机效果.我使用anova函数.我的过程,到目前为止是一个函数调用模型适合lmer()用REML=TRUE(默认选项).然后我打电话anova()给两个模型,其中一个确实包括要测试的随机效果,另一个没有.但是,众所周知,该anova()函数使用ML重新生成模型,但在新版本中,anova()您可以anova()通过设置选项来防止这样做refit=FALSE.为了测试随机效应refit=FALSE,我应该在我的调用中设置anova() or not?(如果我设置refit=FALSEp值往往更低.当我设置时,p值是反保守的refit=FALSE吗？)

方法1:

    mod0_reml <- lmer(x ~ y + z + (1 | w), data=dat)
    mod1_reml <- lmer(x ~ y + z + (y | w), data=dat)
    anova(mod0_reml, mod1_reml)

这将导致anova()使用ML而不是改装模型REML.(该anova()函数的较新版本也将输出有关此信息.)

方法2:

    mod0_reml <- lmer(x ~ y + z + (1 | w), data=dat)
    mod1_reml <- lmer(x ~ y + …

我试图对同一生长季节在2个地点建造的田间试验进行一些统计分析.

在两个站点(Site,水平:HF | NW),实验设计是具有4个(n = 4)块的RCBD(Block每个水平:1 | 2 | 3 | 4 Site).有4种处理--3种不同形式的氮肥和对照(无氮肥)(Treatment水平:AN,U,IU,C).在田间试验期间,有3个不同的时期开始添加肥料,最后收获草.在这个因素下,这些时期被赋予了1 | 2 | 3的水平N_app.

有一系列测量我想测试以下零假设H0:

Treatment (H0)对测量没有影响

我特别感兴趣的两项测量是:草产量和氨排放量.

从这里Dry_tonnes_ha显示的grass yield()开始,一个很好的平衡数据集

可以使用以下代码在R中下载数据:

library(tidyverse)

download.file('https://www.dropbox.com/s/w5ramntwdgpn0e3/HF_NW_grass_yield_data.csv?raw=1', destfile = "HF_NW_grass_yield_data.csv", method = "auto")
raw_data <- read.csv("HF_NW_grass_yield_data.csv", stringsAsFactors = FALSE)

HF_NW_grass <- raw_data %>% mutate_at(vars(Site, N_app, Block, Plot, Treatment), as.factor) %>% 
  mutate(Date = as.Date(Date, format = "%d/%m/%Y"),
         Treatment = factor(Treatment, levels = c("AN", "U", "IU", "C")))

我使用以下方法对此运行ANOVA:

model_1 <- …

我希望能够将命名的模型列表(merMod对象)传递给anova()并在输出中保留模型名称.在使用mclapply()以更高效并行运行一批慢速模型(如glmer)的情况下,这尤其有用.我提出的最好的是在模型列表的去命名版本上使用do.call,但这并不理想,因为我可能有名为(例如)"mod12","mod17"和"mod16"的模型并且这些模型名称在输出中被转换为"MODEL1","MODEL2"和"MODEL3".(在查看单个批次时,这看起来似乎微不足道,但在长时间的建模会话中,有几十个模型,这是一个混乱的确定方法.)

请注意,这与从列表创建和调用线性模型的问题不同,因为我不是要尝试跨列表比较模型对.它比在模型列表上使用lapply更复杂,因为我以非一元的方式使用anova().

这是一个最小的代表:

library(lme4)

formList <- c(mod12 = angle ~ recipe + temp + (1|recipe:replicate),
              mod17 = angle ~ recipe + temperature + (1|recipe:replicate),
              mod16 = angle ~ recipe * temperature + (1|recipe:replicate))
modList <- lapply(formList, FUN=lmer, data=cake)

# Fails because modList is named so it's interpreted as arg-name:arg pairs
do.call(anova, modList)

# Suboptimal because model names aren't preserved
do.call(anova, unname(modList))

# Fails because object isn't merMod (or some other class covered …

我发现statsmodels线性模型的 anova 测试的实现非常有用（http://www.statsmodels.org/dev/generated/statsmodels.stats.anova.anova_lm.html#statsmodels.stats.anova.anova_lm）但我想知道，因为它不存在于库中，如何为逻辑回归部分构建等效版本。

公式：

from statsmodels.formula.api import ols, logit
import statsmodels.api as sm

ols(formula_str, data=data_on_which_to_perform_analysis).fit()
logit(formula_str, data=data_on_which_to_perform_analysis).fit()

sm.stats.anova_lm()

这意味着本质上（通过查看源代码）复制anova_single. 有没有人已经在某个远程存储库中做过这件事？我问是因为实现速度非常快，而且非常深入statsmodels核心库，所以解决它并不容易（至少以我目前的技能水平）

关于如何进行的任何建议？

准备一个小玩具示例：

import pandas as pd
import numpy as np

high, size = 100, 20
df = pd.DataFrame({'perception': np.random.randint(0, high, size),
                   'age': np.random.randint(0, high, size),
                   'outlook': pd.Categorical(np.tile(['positive', 'neutral', 'negative'], size//3+1)[:size]),
                   'smokes': pd.Categorical(np.tile(['lots', 'little', 'not'], size//3+1)[:size]),
                   'outcome': np.random.randint(0, high, size)
                  })
df['age_range'] = pd.Categorical(pd.cut(df.age, range(0, high+5, size//2), right=False,
                             labels=["{0} - {1}".format(i, i + 9) for i in range(0, high, size//2)]))
np.random.shuffle(df['smokes'])

这会给你类似的东西：

In [2]: df.head(10)
Out[2]:
   perception  age   outlook  smokes  outcome age_range
0          13   65  positive  little       22   60 - 69
1 …

我构建了一个Pandas由基因名称索引的数据框（下面的示例），其中包含列和整数的样本名称作为单元格值。我想要做的是对由与样本组对应的列列表定义的行值列表运行 ANOVA ( f_oneway(), from scipy.stats)。然后将这些结果存储在一个新的Pandas数据框中，组名作为列，索引的基因相同。

数据帧的一个例子（它是从 my 中的另一个函数返回的）：

import pandas as pd
counts = {'sample1' : [0, 1, 5, 0, 10],
        'sample2' : [2, 0, 10, 0, 0],
        'sample3' : [0, 0, 0, 1, 0],
        'sample4' : [10, 0, 1, 4, 0]}
data = pd.DataFrame(counts, columns = ['sample1', 'sample2', 'sample3', 'sample4'],
        index = ['gene1', 'gene2', 'gene3', 'gene4', 'gene5'])

组作为参数导入main()，所以在这个函数中我有：

def compare(out_prefix, pops, data):
    import scipy.stats as stats
    sig = pd.DataFrame(index=data.index)

#groups …