使用空间数据从 ddply 函数内部的多元线性回归 (lm) 中提取 p 值

extracting p values from multiple linear regression (lm) inside of a ddply function using spatial data

我有一组空间坐标 (x,y) 数据,其中每个坐标在几年内都有一个响应变量。以下代码生成类似的数据框:

df <- data.frame(
  id = rep(1:2, 2),
  x = rep(c(25, 30),10),
  y = rep(c(100, 200), 10),
  year = rep(1980:1989, 2),
  response = rnorm(20)
)

结果数据框:

head(df)

 id  x   y year   response
1  1 25 100 1980  0.1707431
2  2 30 200 1981  1.3562263
3  1 25 100 1982 -0.4590506
4  2 30 200 1983  1.3238410
5  1 25 100 1984  1.7765772
6  2 30 200 1985 -0.6258069

我想 运行 随时间对每个单元格进行线性回归以获得响应变量的变化。为此,我使用 plyr 和 ddply 命令:

require(plyr)    
lm.df <- ddply(df, .(id), function(z)coef(lm(response ~ year, data = z)))

然后我重新组合空间数据(这里的重新组合示例有点简单,但适用于示例。):

points <- data.frame(
  id = c(1,2),
  x = c(25,30),
  y = c(100,200)
)

lm.stack <-  merge(points, lm.df, by="id")
colnames(lm.stack) <- c("ID", "x", "y", "intercept", "slope")

print(lm.stack)
ID  x   y intercept       slope
1  1 25 100  257.7291 -0.12985632
2  2 30 200  173.3676 -0.08708068

效果很好。但我想要的是能够提取每个 lm 模型的 adj r 平方值,以便能够通过坐标单元表示显着趋势,最终目标是投影按显着性值着色的单元格图。这是我需要帮助的地方,我非常感激。谢谢!

如果您将模型保存在列表中,那么您可以使用 lapply 或任何其他列表遍历函数轻松地从中提取您想要的任何内容。

lms <- dlply(df, .(id), function(x) lm(response ~ year, data=x))
rsq <- lapply(lms, function(x) summary(x)$adj.r.squared)  # get the adjusted r2

您自己制作列表的另一种方法是 lmList 来自 nlme

library(nlme)
lms2 <- lmList(response ~ year | id, data=df)
summary(lms2)$adj.r.squared

我建议组合使用包 dplyrbroom。 这种方法将 return 您需要了解的有关模型的所有信息(有关系数的信息和有关模型本身的信息)作为数据框:

set.seed(25)

df <- data.frame(
  id = rep(1:2, 2),
  x = rep(c(25, 30),10),
  y = rep(c(100, 200), 10),
  year = rep(1980:1989, 2),
  response = rnorm(20)
)

df

#    id  x   y year    response
# 1   1 25 100 1980 -0.21183360
# 2   2 30 200 1981 -1.04159113
# 3   1 25 100 1982 -1.15330756
# 4   2 30 200 1983  0.32153150
# 5   1 25 100 1984 -1.50012988
# 6   2 30 200 1985 -0.44553326
# 7   1 25 100 1986  1.73404543
# 8   2 30 200 1987  0.51129562
# 9   1 25 100 1988  0.09964504
# 10  2 30 200 1989 -0.05789111
# 11  1 25 100 1980 -1.74278763
# 12  2 30 200 1981 -1.32495298
# 13  1 25 100 1982 -0.54793388
# 14  2 30 200 1983 -1.45638428
# 15  1 25 100 1984  0.08268682
# 16  2 30 200 1985  0.92757895
# 17  1 25 100 1986 -0.71676933
# 18  2 30 200 1987  0.96239968
# 19  1 25 100 1988  1.54588458
# 20  2 30 200 1989 -1.00976361



library(dplyr)
library(broom)


df %>% 
  group_by(id) %>%
  do({model = lm(response~year, data=.)    # create your model
      data.frame(tidy(model),              # get coefficient info
                 glance(model))})          # get model info


#   id        term     estimate    std.error statistic    p.value r.squared adj.r.squared     sigma statistic.1  p.value.1 df    logLik      AIC      BIC deviance df.residual
# 1  1 (Intercept) -492.2144842 213.19252113 -2.308779 0.04978386 0.3996362    0.32459069 0.9611139    5.325253 0.04987162  2 -12.67704 31.35409 32.26184 7.389919           8
# 2  1        year    0.2479705   0.10745580  2.307651 0.04987162 0.3996362    0.32459069 0.9611139    5.325253 0.04987162  2 -12.67704 31.35409 32.26184 7.389919           8
# 3  2 (Intercept) -258.6253012 196.88284243 -1.313600 0.22539989 0.1771294    0.07427055 0.8871395    1.722063 0.22582607  2 -11.87614 29.75227 30.66003 6.296132           8
# 4  2        year    0.1301582   0.09918521  1.312274 0.22582607 0.1771294    0.07427055 0.8871395    1.722063 0.22582607  2 -11.87614 29.75227 30.66003 6.296132           8

如果您真的想要最终数据框中的 x 和 y 列,则可以使用 group_by(id,x,y)。例如,如果您想要与您提供的输出类似的输出,您可以这样做:

library(dplyr)
library(broom)
library(tidyr)

dd = 
    df %>% 
      group_by(id, x, y) %>%
      do({model = lm(response~year, data=.)    # create your model
      data.frame(tidy(model),              # get coefficient info
                 glance(model))})          # get model info

然后select您需要的信息:

 dd %>%
  select(id, x, y, term, estimate, adj.r.squared)

#   id  x   y        term     estimate adj.r.squared
# 1  1 25 100 (Intercept) -492.2144842    0.32459069
# 2  1 25 100        year    0.2479705    0.32459069
# 3  2 30 200 (Intercept) -258.6253012    0.07427055
# 4  2 30 200        year    0.1301582    0.07427055

截距一行,斜率一行。

或者,甚至将该数据框重塑为:

dd %>%
  select(id, x, y, term, estimate, adj.r.squared) %>%
  spread(term, estimate)

#   id  x   y adj.r.squared (Intercept)      year
# 1  1 25 100    0.32459069   -492.2145 0.2479705
# 2  2 30 200    0.07427055   -258.6253 0.1301582

yearslope(变量系数year

以类似的方式,您还可以使用包 purrr 创建一个新的数据框,其中包含带有相应信息的列表列:

set.seed(25)

df <- data.frame(
  id = rep(1:2, 2),
  x = rep(c(25, 30),10),
  y = rep(c(100, 200), 10),
  year = rep(1980:1989, 2),
  response = rnorm(20)
)

library(dplyr)
library(broom)
library(tidyr)
library(purrr)

dd = df %>%
  group_by(id, x, y) %>%                                  # for each combination of those variables
  nest() %>%                                              # nest data (rest of columns)
  mutate(Model = map(data, ~lm(response~year, data=.)),   # use nested data to build model of interest
         Coeff_Info = map(Model, tidy),                   # get coefficient info
         Model_Info = map(Model, glance)) %>%             # get model info
  ungroup()                                               # forget the grouping

# check how new dataset looks like
dd

# # A tibble: 2 x 7
#        id     x     y              data    Model           Coeff_Info            Model_Info
#     <int> <dbl> <dbl>            <list>   <list>               <list>                <list>
#   1     1    25   100 <tibble [10 x 2]> <S3: lm> <data.frame [2 x 5]> <data.frame [1 x 11]>
#   2     2    30   200 <tibble [10 x 2]> <S3: lm> <data.frame [2 x 5]> <data.frame [1 x 11]>

您仍然可以访问任何您想要的元素,但请记住,您的一些列现在有列表元素:

# get coefficient info for both models
dd$Coeff_Info

# [[1]]
#          term     estimate   std.error statistic    p.value
# 1 (Intercept) -492.2144842 213.1925211 -2.308779 0.04978386
# 2        year    0.2479705   0.1074558  2.307651 0.04987162
# 
# [[2]]
#          term     estimate    std.error statistic   p.value
# 1 (Intercept) -258.6253012 196.88284243 -1.313600 0.2253999
# 2        year    0.1301582   0.09918521  1.312274 0.2258261


# get the r squared value for 1st model
dd %>%                    # from new dataset
  filter(id == 1) %>%     # keep all info / rows of 1st model
  pull(Model_Info) %>%    # get model info
  map_dbl("r.squared")    # show r squared

# [1] 0.3996362

或者,或者,

dd %>%                        # from new dataset
  unnest(id, Model_Info) %>%  # umnest id and model info
  filter(id == 1) %>%         # keep row of 1st model
  pull(r.squared)             # show the r squared

# [1] 0.3996362