Extracting final p-value statistic from an lm lapply loop with multiple models

user_15431348 0
Viewed 293

I have the following code that automatically performs lm between my independent variable (Kpl) and all my other dependent variables (Y1, Y2, ...., Yi):

linear_summary <- lapply(testdata[,-1], function(x) summary(lm(Kpl ~ x)))

The output for this is


Call:
lm(formula = Kpl ~ x)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.37567 -0.52392  0.04236  0.67444  0.81316 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.7282     0.3456   5.001 0.000402 ***
x            -0.1550     0.2712  -0.571 0.579196    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.772 on 11 degrees of freedom
Multiple R-squared:  0.02883,   Adjusted R-squared:  -0.05946 
F-statistic: 0.3265 on 1 and 11 DF,  p-value: 0.5792


$Y2

Call:
lm(formula = Kpl ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.2472 -0.4236 -0.2057  0.7140  1.0348 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   0.6900     0.9010   0.766    0.460
x             0.8832     0.8767   1.007    0.335

Residual standard error: 0.7495 on 11 degrees of freedom
Multiple R-squared:  0.08447,   Adjusted R-squared:  0.001238 
F-statistic: 1.015 on 1 and 11 DF,  p-value: 0.3354

Etc. (I have truncated it for just the first 2 correlations)

I wanted to extract the final p-value for the whole model for each of the instances (0.5792 and 0.3354 in these two cases). Ideally this would come in some sort of table form with the associated correlation variable, i.e. Y1=0.5792 Y2=0.3354.

Most of the info I can find either seem to only work for a single correlation (as opposed to an sapply with multiple correlations) or I do not seem to get it to work, which could be a problem with my original code.

Any suggestions for a person just starting with R on how to solve this?

Edit: The data looks something like this

|    X     |     Y1      |     Y2      |     Y3      |     Y4      |
| -------- | ------------|-------------|-------------|-------------|
| 0.33767  | 2.33063062  | 1.013212308 | 1.277996888 | 1.373238355 |
| 0.33767  | 0.095967324 | 0.508830529 | 0.789257027 | 0.815877121 |
| 1.010474 | 2.344657045 | 0.842490752 | 1.240582283 | 1.262360905 |
| 1.010474 | 0.08135992  | 0.912535398 | 0.384427466 | 0.409817599 |
| 1.183276 | 0.135626937 | 0.967877981 | 0.505801442 | 0.576288093 |
| 1.536974 | 1.507146148 | 1.428839993 | 1.316569449 | 1.392022619 |
| 1.536974 | 1.255210981 | 1.191822955 | 1.395769591 | 1.41903939  |
| 2.017965 | 1.410299711 | 1.121560244 | 1.369835675 | 1.385143026 |
| 2.017965 | 1.032587109 | 1.372235121 | 1.390878783 | 1.42741762  |
| 2.3436   | 1.275999998 | 0.930400789 | 1.19877482  | 1.217540034 |
| 2.3436   | 1.250513383 | 1.063880146 | 1.206719195 | 1.23325973  |
| 2.387598 | 0.182866909 | 0.89588293  | 0.416923749 | 0.45364797  |
| 2.387598 | 0.097133916 | 0.750430855 | 0.506463633 | 0.03434754  |

These are the actual values that I used to get the correlations above

sapply lm linear-regression r
2 Answers

I think the p-value is not stored, you need to calculate it from the fstatistics, maybe something like this:

set.seed(111)
testdata = data.frame(Kpl = rnorm(100), Y1 = rnorm(100),
                      Y2 = rnorm(100), Y3 = rnorm(100))

IV = colnames(testdata)[-1]
DV = "Kpl"

linear_summary <- lapply(IV,function(x){
         summary(lm(reformulate(response=DV,termlabels=x),data=testdata))
                         })

names(linear_summary) = IV

tab = lapply(IV,function(x){
  p = with(
       linear_summary[[x]],
       pf(fstatistic[1],fstatistic[2],fstatistic[3],lower.tail=FALSE)
          )
  data.frame(IV = x, p = p)
})

do.call(rbind,tab)

       IV         p
value  Y1 0.5757187
value1 Y2 0.4922582
value2 Y3 0.4009439

Check for example first summary:

linear_summary[[1]]

Call:
lm(formula = reformulate(response = DV, termlabels = x), data = testdata)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.94515 -0.73325  0.05448  0.57901  2.76026 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.01382    0.10747  -0.129    0.898
Y1          -0.05950    0.10597  -0.562    0.576

Residual standard error: 1.075 on 98 degrees of freedom
Multiple R-squared:  0.003207,  Adjusted R-squared:  -0.006964 
F-statistic: 0.3153 on 1 and 98 DF,  p-value: 0.5757
user_122584590

Ok I edited my code in the following way:

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

df %>%    # Solution 1
  pivot_longer(-X) %>%
  group_split(name) %>%
  set_names(nm = map(., ~ first(.x$name))) %>%
  map(~ tidy(lm(X ~ value, data = .))) %>%
  bind_rows(.id = "var") %>%
  filter(term == "value")

# A tibble: 4 x 6
  var   term  estimate std.error statistic p.value
  <chr> <chr>    <dbl>     <dbl>     <dbl>   <dbl>
1 Y1    value  -0.155      0.271   -0.571    0.579
2 Y2    value   0.883      0.877    1.01     0.335
3 Y3    value   0.0341     0.552    0.0618   0.952
4 Y4    value  -0.158      0.469   -0.337    0.743

Or you can use this:

df %>%    # Solution 2
  pivot_longer(Y1:Y4) %>%
  group_by(name) %>%
  arrange(.by_group = TRUE) %>% 
  nest() %>%
  mutate(models = map(data, ~ lm(X ~ value, data = .)),
         glance = map(models, glance)) %>%
  unnest(glance)

# A tibble: 4 x 15
# Groups:   name [4]
  name  data    models r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
  <chr> <list>  <list>     <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl>
1 Y1    <tibbl~ <lm>    0.0288        -0.0595  0.772   0.327     0.579     1  -14.0  34.0  35.7
2 Y2    <tibbl~ <lm>    0.0845         0.00124 0.750   1.01      0.335     1  -13.6  33.2  34.9
3 Y3    <tibbl~ <lm>    0.000348      -0.0905  0.783   0.00382   0.952     1  -14.2  34.4  36.1
4 Y4    <tibbl~ <lm>    0.0102        -0.0798  0.779   0.113     0.743     1  -14.1  34.2  35.9
# ... with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

I know you already got your answer but here I presented 2 other solutions. Thought it may be ok to learn alternative ways of dealing with a problem and thank you for your question, it was very good.

user_143145200
Related