多次拟合回归并收集汇总统计信息答案

【问题标题】：Fitting regression multiple times and gather summary statistics多次拟合回归并收集汇总统计信息
【发布时间】：2019-06-06 00:02:38
【问题描述】：

我有一个如下所示的数据框：

W01           0.750000     0.916667     0.642857      1.000000      0.619565   
W02           0.880000     0.944444     0.500000      0.991228      0.675439   
W03           0.729167     0.900000     0.444444      1.000000      0.611111   
W04           0.809524     0.869565     0.500000      1.000000      0.709091   
W05           0.625000     0.925926     0.653846      1.000000      0.589286   

Variation  1_941119_A/G  1_942335_C/G  1_942451_T/C  1_942934_G/C  \
W01            0.967391      0.965909             1      0.130435   
W02            0.929825      0.937500             1      0.184211   
W03            0.925926      0.880000             1      0.138889   
W04            0.918182      0.907407             1      0.200000   
W05            0.901786      0.858491             1      0.178571   

Variation  1_944296_G/A    ...     X_155545046_C/T  X_155774775_G/T  \
W01            0.978261    ...            0.652174         0.641304   
W02            0.938596    ...            0.728070         0.736842   
W03            0.944444    ...            0.675926         0.685185   
W04            0.927273    ...            0.800000         0.690909   
W05            0.901786    ...            0.794643         0.705357   

Variation  Y_5100327_G/T  Y_5100614_T/G  Y_12786160_G/A  Y_12914512_C/A  \
W01             0.807692       0.800000        0.730769        0.807692   
W02             0.655172       0.653846        0.551724        0.666667   
W03             0.880000       0.909091        0.833333        0.916667   
W04             0.666667       0.642857        0.580645        0.678571   
W05             0.730769       0.720000        0.692308        0.720000   

Variation  Y_13470103_G/A  Y_19705901_A/G  Y_20587967_A/C  mean_age  
W01              0.807692        0.666667        0.333333      56.3  
W02              0.678571        0.520000        0.250000      66.3  
W03              0.916667        0.764706        0.291667      69.7  
W04              0.666667        0.560000        0.322581      71.6  
W05              0.703704        0.600000        0.346154      72.5  

[5 rows x 67000 columns]

我想为每列拟合一个简单的最小二乘线性回归和 Thiel-Sen 线性回归作为自变量，平均年龄作为响应变量，并收集汇总统计数据，包括 slope、intercept、@ 987654327@、p value 和 std err 用于每个拟合，最好将输出收集为数据农场！

到目前为止，我一直在对“df”进行切片并分别对每一列进行回归分析：

from scipy import stats
import time

# Start timer
start_time = time.time()

# Select only 'Variation of interest' and 'mean_age' columns
r1 = tdf [['1_944296_G/A', 'mean_age']]

# Use scipy lingress function to perform linear regression
slope, intercept, r_value, p_value, std_err = stats.linregress(tdf['mean_age'], \
    tdf['1_69270_A/G'])

print('The p-value between the 2 variables is measured as ' + str(p_value) + '\n')
print('Least squares linear model coefficients, intercept = ' + str(intercept) + \
  '. Slope = ' + str(slope)+'\n')

# Create regression line
regressLine = intercept + tdf['mean_age']*slope

# Regression using Theil-Sen with 95% confidence intervals 
res = stats.theilslopes(tdf['1_69270_A/G'], tdf['mean_age'], 0.95)

print('Thiel-Sen linear model coefficients, intercept = ' + str(res[1]) + '. Slope = ' + \
  str(res[0]) +'\n')

# Scatter plot the temperature
plt.clf()
plt.scatter(tdf['mean_age'], tdf['1_69270_A/G'], s = 3, label = 'Allele frequency')

# Add least squares regression line
plt.plot(tdf['mean_age'], regressLine, label = 'Least squares regression line'); 

# Add Theil-Sen regression line
plt.plot(tdf['mean_age'], res[1] + res[0] * tdf['mean_age'], 'r-', label = 'Theil-Sen regression line')

# Add Theil-Sen confidence intervals
plt.plot(tdf['mean_age'], res[1] + res[2] * tdf['mean_age'], 'r--', label = 'Theil-Sen 95% confidence interval')
plt.plot(tdf['mean_age'], res[1] + res[3] * tdf['mean_age'], 'r--')

# Add legend, axis limits and save to png
plt.legend(loc = 'upper left')
#plt.ylim(7,14); plt.xlim(1755, 2016)
plt.xlabel('Year'); plt.ylabel('Temperature (C)')
plt.savefig('pythonRegress.png')

# End timer
end_time = time.time()
print('Elapsed time = ' + str(end_time - start_time) + ' seconds')

我想知道如何在每个列的迭代循环中执行此分析，并将最终结果收集到一个综合数据框中。

我见过[this](Looping regression and obtaining summary statistics in matrix form"循环回归，获取矩阵形式的汇总统计 ")！但不是我期望的输出。Python 或 R 中的任何解决方案都值得赞赏！

【问题讨论】：

请设置一个minimal，可运行，reproducible example。如果您询问摘要统计信息，为什么要包括绘图线？请将问题集中在 minimal 示例上。

标签： python r pandas linear-regression

【解决方案1】：

我认为您会发现本指南很有用：Running a model on separate groups。

让我们生成一些与您的数据相似的示例数据，其中包含两个变量和平均年龄的值。我们还需要几个包：

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

set.seed(1001)
data1 <- data.frame(mean_age = sample(40:80, 50, replace = TRUE), 
                    snp01 = rnorm(50), 
                    snp02 = rnorm(50))

第一步是使用gather 将“宽”格式转换为“长”格式，以便变体名称在一个列中，而值在另一列中。然后我们可以通过变体名称nest。

data1 %>% 
  gather(snp, value, -mean_age) %>% 
  nest(-snp)

这将创建一个 tibble（一个特殊的数据框），其中第二列 data 是一个“列表列” - 它包含平均年龄和该行中变量的值：

# A tibble: 2 x 2
  snp   data             
  <chr> <list>           
1 snp01 <tibble [50 x 2]>
2 snp02 <tibble [50 x 2]>

现在我们使用purrr::map 为每一行创建带有线性模型的第三列：

data1 %>% 
  gather(snp, value, -mean_age) %>% 
  nest(-snp) %>% 
  mutate(model = map(data, ~lm(mean_age ~ value, data = .)))

结果：

# A tibble: 2 x 3
  snp   data              model 
  <chr> <list>            <list>
1 snp01 <tibble [50 x 2]> <lm>  
2 snp02 <tibble [50 x 2]> <lm>

最后一步是根据需要汇总模型，然后unnest 数据结构。我正在使用broom::glance()。完整流程：

data1 %>% 
  gather(snp, value, -mean_age) %>% 
  nest(-snp) %>% 
  mutate(model = map(data, ~lm(mean_age ~ value, data = .)), 
         summary = map(model, glance)) %>% 
  select(-data, -model) %>% 
  unnest(summary)

结果：

# A tibble: 2 x 12
  snp   r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC deviance df.residual
  <chr>     <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>  <dbl> <dbl> <dbl>    <dbl>       <int>
1 snp01   0.00732      -0.0134   12.0     0.354   0.555     2  -194.  394.  400.    6901.          48
2 snp02   0.0108       -0.00981  12.0     0.524   0.473     2  -194.  394.  400.    6877.          48

【讨论】：

有什么办法也可以得到截距值吗？
你可以使用broom::tidy()而不是glance来获取截取值。
是否可以在汇总输出中也包含回归斜率？
broom::tidy 返回估计值，即预测变量的斜率。
尼尔，我想知道你能否帮助我解决另一个问题。我想将 robust::lmrob() 与 MM-estimator 一起使用，但一目了然或整洁不支持。您对我如何收集与上面的 sn-p 一样简洁的汇总统计数据有什么建议吗？

【解决方案2】：

我不知道您的数据和分析的确切细节和复杂性，但这是我会采取的方法。

data <- data.frame(mean_age=rnorm(5),
                   Column_1=rnorm(5),
                   Column_2=rnorm(5),
                   Column_3=rnorm(5),
                   Column_4=rnorm(5),
                   Column_5=rnorm(5)
                   )
data


looped <- list()

for(each_col in names(data)[-1]){
    looped[[each_col]] <- lm(get(each_col) ~ mean_age, data)

}

looped

【讨论】：