尝试 运行 数据帧中字段的方差分析 [R]

Trying to run an ANOVA for fields in a dataframe [R]

到目前为止,我的代码如下所示:

Points = readOGR(dsn = "./Data/filename.shp",layer = "layername",stringsAsFactors = FALSE) 
Points$LDI = extract(LDI, Points) 
LDI = raster("./Data/filename2.tif")
Points$LDI = extract(LDI, Points)
PointsDF = Points@data

for(i in PointsDF) { 
    Mod1 = lm(LDI ~ i, data = PointsDF) 
    Mod2 = lm(LDI ~ 1, data = PointsDF) 
    anova(Mod1, Mod2)
}

最后一部分是我知道我做错的地方。我想 运行 数据框中每个数字字段的方差分析。

问题是您没有告诉循环它正在迭代什么,没有在方差分析调用中定义公式对象,也没有创建对象来存储结果。

在此示例中,"ij" 变量分配给列表对象并存储方差分析模型,"y" 定义为指示模型 left-hand 侧的变量。列表对象 "anova.results" 正在存储每个模型。循环定义中的索引使用 "which" 来分配包含 "y" 的列,因此将其从迭代器中删除。我正在使用 R "iris" 数据集作为示例。

data(iris)
iris <- iris[,-5]

y = "Sepal.Length"    
anova.results <- list()
ij=0
  for(i in names(iris)[-which(names(iris) == y)]) {
    ij = ij+1
    Mod = lm(stats::as.formula(paste(y, i, sep = "~")), data = iris) 
    anova.results[[ij]] <- anova(Mod, Mod)
  }
anova.results

你很接近。一种自然的方法是遍历字段名称。虽然有很多方法可以做到这一点,但 lapply 可能是最惯用的,因为 (a) 它使用字段名称(而不是字段索引,这可能很危险)并且 (b) 不需要 pre-allocating 输出的任何结构。诀窍是将字段名称转换为公式。同样,有很多方法可以做到这一点,但直接的方法是 assemble 将公式作为字符串。

这里以工作代码为例。它生成 anova 个对象的列表。

#
# Create some random data.
#
n <- 20
set.seed(17)
X <- data.frame(Y=rnorm(n), X1=runif(n), X2=1:n, X3=rexp(n))
#
# Loop over the regressors.
# (The base model can be precomputed.)
#
mod.0 <- lm(Y ~ 1, X)
models <- lapply(setdiff(names(X), "Y"), function(s) {
  mod.1 <- lm(as.formula(paste("Y ~", s)), X)
  anova(mod.0, mod.1)
})
print(models)

这是输出,显示了三个 anova 结果的列表。

[[1]]
Analysis of Variance Table

Model 1: Y ~ 1
Model 2: Y ~ X1
  Res.Df     RSS Df Sum of Sq     F Pr(>F)
1     19 10.1157                          
2     18  9.6719  1   0.44385 0.826 0.3754

[[2]]
Analysis of Variance Table

Model 1: Y ~ 1
Model 2: Y ~ X2
  Res.Df     RSS Df Sum of Sq      F  Pr(>F)  
1     19 10.1157                              
2     18  8.1768  1     1.939 4.2684 0.05353 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

[[3]]
Analysis of Variance Table

Model 1: Y ~ 1
Model 2: Y ~ X3
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     19 10.116                           
2     18 10.081  1  0.034925 0.0624 0.8056

作为使用您制作的内容的另一个示例,这里是 sapply 用于打印他们的 p-values:

sapply(models, function(m) m[["Pr(>F)"]][2])

[1] 0.37542968 0.05352883 0.80562894