尝试：

results <- list()
for(i in names(ados.sim[,2:4])){  
  results[[i]] <- kruskal.test(formula(paste(i, "~ muttype")), data = ados.sim)
}

这也将结果保存在列表中，避免在每次迭代中将结果覆盖为ados.mtp，我认为这不是您想要做的。
请注意以下事项：

for(i in names(ados.sim[,2:4])){  
   print(i)
}
[1] "ados.tsc"
[1] "ados.sa"
[1] "ados.rrb"

也就是说，i已经给出了列的名称。代码中的问题是，您试图将其用作整数来进行子集设置，这将结果变为NA。

for(i in names(ados.sim[,2:4])){  
   print(paste((colnames(ados.sim)[i]), "~ muttype"))
}
[1] "NA ~ muttype"
[1] "NA ~ muttype"
[1] "NA ~ muttype"

仅供参考，所有这些也可以通过以下两种我经常喜欢的方式来完成，因为它使后续分析稍微容易一些：
首先，将所有测试对象存储在 Dataframe 中：

library(tidyr)
df <- ados.sim %>% gather(key, value, -muttype) %>% 
      group_by(key) %>% 
      do(test = kruskal.test(x= .$value, g = .$muttype))

然后，您可以对dataframe进行子集设置以获得测试结果：

df[df$key == "ados.rrb",]$test
[[1]]

    Kruskal-Wallis rank sum test

data:  .$value and .$muttype
Kruskal-Wallis chi-squared = 2.2205, df = 2, p-value = 0.3295

或者，直接在dataframe中获取所有结果，而不存储测试对象：

library(broom)
df2 <- ados.sim %>% gather(key, value, -muttype) %>% 
       group_by(key) %>% 
       do(tidy(kruskal.test(x= .$value, g = .$muttype)))
df2
# A tibble: 3 x 5
# Groups:   key [3]
       key statistic   p.value parameter                       method
     <chr>     <dbl>     <dbl>     <int>                       <fctr>
1 ados.rrb 2.2205031 0.3294761         2 Kruskal-Wallis rank sum test
2  ados.sa 0.1319554 0.9361517         2 Kruskal-Wallis rank sum test
3 ados.tsc 0.3618102 0.8345146         2 Kruskal-Wallis rank sum test

如果可以使用外部库：

library(matrixTests)

col_kruskalwallis(ados.sim[,-1], ados.sim[,1])

         obs.tot obs.groups df statistic    pvalue
ados.tsc     240          3  2 0.3618102 0.8345146
ados.sa      240          3  2 0.1319554 0.9361517
ados.rrb     240          3  2 2.2205031 0.3294761