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

Question

到目前为止，我的代码如下所示：

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)
}

最后一部分是我知道我做错了一切的地方。我想在数据帧中的每个数值字段上运行方差分析。

Answer 1

你很接近。 一种自然的方法是遍历字段名称。 尽管有很多方法可以做到这一点，但lapply可能是最惯用的，因为 （a） 它使用字段名称（而不是字段索引，这可能是危险的）和 （b） 不需要为输出预先分配任何结构。 诀窍是将字段名称转换为公式。 同样，有很多方法可以做到这一点，但直接的方法是将公式组装为字符串。

下面是工作代码作为示例。 它生成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

作为处理您生成的内容的另一个示例，以下是用于打印其 p 值sapply：

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

[1] 0.37542968 0.05352883 0.80562894

Answer 2

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

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

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