使用R语言的`modifiedmk`包中的统计量来创建汇总表格

Question

使用R语言的`modifiedmk`包中的统计量来创建汇总表格

5

我正在尝试在R中运行来自modifiedmk包的函数。

install.packages('modifiedmk')
library(modifiedmk)

我有一个数据框 data，是通过以下方式生成的：

Station <- c('APT','APT', 'APT','APT', 'APT', 'APT', 'APT','APT', 'APT','APT','APT','APT',
              'AF','AF', 'AF','AF','AF','AF','AF','AF','AF',
             'EL', 'EL', 'EL', 'EL', 'EL', 'EL', 'EL', 'EL', 'EL', 'EL', 'EL', 'EL', 'EL', 'EL', 'EL',
             'GFS', 'GFS', 'GFS', 'GFS', 'GFS', 'GFS', 'GFS', 'GFS', 'GFS', 'GFS', 'GFS', 'GFS', 'GFS', 'GFS', 'GFS', 'GFS'
              )
Rainfall <- c(375.3, 263.3, 399.2, 242.6, 847.6, 276.5, 712.8, 366.3, 188.6, 478.4, 539, 682.5,
            520.7, 1337.8, 524, 908.4,748.5,411.8, 772.4,978.5,983,
            732.4, 788.6, 567.1, 576, 931.6, 727.2, 1079.3, 902.8,493.4,  630.7, 784.1,660.2, 531.3, 487.1,798.4,
            1064.1,  590.3, 1011.2, 1037.1,  1398.4, 1153.6,994.1,  1100.2,743.7,637.4, 792.2, 891.9,880.9, 670, 920.2,681.4)
Year <- c('1957','1958','1959','1960','1961','1962','1963','1964','1965','1966','1967','1968',
                  '1960','1961','1962','1963','1964','1965','1966','1967','1968',
                  '1957','1958','1959','1960','1961','1962','1963','1964','1965','1966','1967','1968','1969','1970','1971',
                  '1964','1965','1966','1967','1968','1969','1970','1971','1972','1973','1974','1975','1976','1977','1978','1979')
length(Year)
data<-data.frame(Year, Station, Rainfall)

我有一个数据框，其中四个降雨数据的站点作为行。我想在每个数据的 Station 上应用来自 modifiedmk 软件包的 mmky1lag 方法，并在R中生成具有两列的摘要表格：

在p <0.05的显着趋势的站点的百分比
Sen斜率的平均值

例如，我可以使用 mmky1lag（as.vector（data $ Rainfall））对所有降雨数据运行 mmky1lag 方法，它会产生

> mmky1lag(as.vector(data$Rainfall))
Corrected Zc  new P-value         N/N*   Original Z  old P.value 
3.332353e+00 8.611480e-04 1.297360e+00 3.795608e+00 1.472822e-04 
         Tau  Sen's slope old.variance new.variance 
3.634992e-01 9.092857e+00 1.605933e+04 2.083474e+04

"我对这些输出中的两个感兴趣："

"第一列："

# Get percent of stations with significant trends where p < 0.05
mmky1lag(as.vector(data$Rainfall))[2] < 0.05

"和第二列："

# Make another column that is the mean Sen's slope
mmky1lag(as.vector(data$Rainfall))[7]

然而，我如何将此方法应用于每个单独的Station的data中？在Python中，我会按Station进行分组，然后应用该方法。但我不确定如何在R中实现这一点。然后，在按站点分组之后，我想要一个包含上述两列信息的汇总表格。

- JAG2024

4个回答

2

你可以尝试像这样在基本R中完成一些操作。首先，你可以将数据作为列表，每个元素都是一个“Station”。

data_list <- split(data,data$Station)

你可以使用lapply()，引自文档:

lapply返回一个与X长度相同的列表，其中每个元素是将FUN应用于X对应元素的结果。

library(modifiedmk)
stat_list <- lapply(data_list, function(x) mmky1lag(x$Rainfall))

现在，你可以将其作为data.frame，然后计算你需要的内容。你可以使用do.call()将rbind()应用到列表中，并将其放入data.frame()中。通常我更喜欢使用列名而不是它们的索引进行操作，但这是主观的。从文档rbind()： “按列或行将向量、矩阵或数据框序列组合起来。这些都是具有其他R类方法的通用函数。” 从文档do.call() ： “从名称或函数和参数列表构造并执行函数调用。”

stat_df <- data.frame(do.call(rbind, stat_list))

现在您可以轻松计算所需内容：

# percentage of the < 0.05 p-values
# here you calculate the number of row of the subset of interest of the
# df / number of row of the dataset.
nrow(stat_df[stat_df$new.P.value < 0.05,])/nrow(stat_df)*100
[1] 0

# Or if you want a prettier result printed:
library(formattable)
percent(nrow(stat_df[stat_df$new.P.value < 0.05,])/nrow(stat_df))
[1] 0.00%

# the mean of Sen.s.slope
mean(stat_df$Sen.s.slope)
[1] 4.45125

此外，我不明白您希望的期望输出方式，它写着“Column1”和“Column2”。如果您定义清楚，就可以得到更符合您要求的结果。

- s__

1

这个符合要求吗？百分比将为零，因为所有的p值都大于5%。你需要在循环中添加< 0.05才能在数据框中获得真/假值。

results <- data.frame(matrix(NA, 4, 3))
colnames(results) <- c('station', 'p-val', 'Sen-slope')
for(ii in seq_along(unique(Station))){
  i <- unique(Station)[ii]
  results[ii, 1] <- i
  results[ii, 2] <- mmky1lag(as.vector(data$Rainfall[data$Station %in% i]))[2]
  results[ii, 3] <- mmky1lag(as.vector(data$Rainfall[data$Station %in% i]))[7]
}

> results
  station     p-val Sen-slope
1     APT 0.2005025  17.32083
2      AF 0.2140574  29.73750
3      EL 0.4561249  -9.60000
4     GFS 0.1854354 -19.65333

- tester

0

如果您使用tidyverse，那么像pandas一样的语法就很容易实现。

# Importing tidyverse
library(tidyverse)

# Calculating grouped values
data %>%
  group_by(Station) %>%
  summarise('p-value' = mmky1lag(Rainfall)[2]<0.05, "Sen's slope" = mmky1lag(Rainfall)[7])

# Output
Station p.value Sen.s.slope
AF      FALSE    29.73750
APT     FALSE    17.32083
EL      FALSE    -9.60000
GFS     FALSE   -19.65333

- Esben Eickhardt

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ben · Accepted Answer

如果你想按组（在这种情况下是“station”）对数据框应用“mmky1lag”函数，有多种方法可供考虑。

首先，你可以使用“aggregate”：

library(modifiedmk)

mktests <- aggregate(Rainfall ~ Station, data = data, FUN = mmky1lag)

这将使用Station组的Rainfall度量来计算一个公式。您的所有结果将以矩阵形式返回，MK测试参数将在单列中。

另一种方法可能是使用data.table包。

library(data.table)

mktests <- as.data.table(data)[, as.list(mmky1lag(Rainfall)), by = Station]

这将把mmky1lag的结果放入列表中，并转换为数据表。选项by允许您按Station执行此操作。第三种方法可能涉及使用dplyr包。

library(dplyr)

mktests <- data %>%
  group_by(Station) %>%
  group_map(~mmky1lag(.x$Rainfall)) %>%
  setNames(unique(sort(data$Station))) %>%
  bind_rows(.id = "Station")

这里使用 group_by 按照 Station 进行分组，然后使用 group_map 对分组元素应用 mmky1lag 函数。使用 setNames 将 Station 的值添加回结果中，再使用 bind_rows 将结果列表转换为数据框。

使用 data.table 解决方案得到的结果如下（其他方法的结果类似）：

R> mktests
   Station Corrected Zc new P-value      N/N* Original Z old P.value        Tau Sen's slope old.variance new.variance
1:     APT    1.2801214   0.2005025 0.4849366  0.8914431   0.3726915  0.2121212    17.32083     212.6667    103.12986
2:      AF    1.2424858   0.2140574 0.5703144  0.9383149   0.3480826  0.2777778    29.73750      92.0000     52.46892
3:      EL   -0.7452428   0.4561249 1.1288325 -0.7917947   0.4284804 -0.1619048    -9.60000     408.3333    460.93994
4:     GFS   -1.3242038   0.1854354 1.4160741 -1.5757881   0.1150746 -0.3000000   -19.65333     493.3333    698.59657

如果你想要计算 p 值小于 0.05 的 Station 百分比，你可以这样做：

sum(mktests$`new P-value` < .05) / nrow(mktests)

在这种情况下，由于基于新的P值，它们中没有一个是显著的，因此为零。

可以计算出Sen的斜率的平均值：

mean(mktests$`Sen's slope`)
4.45125

我不确定您是否预期使用示例数据会得到不同的结果（因为您建议结果将放入2列）。如果这是您想要的，请告诉我。