今天被我发现线性回归之后的众多系数都可以秒秒钟提取的呀，summary.lm让科研更方便，一定要记下来。
**************************************************************************************************************
现在有一批气象站点的数据，某研究想看一下站点的温度与高程和经纬度的关系，其实重点是温度随高程的变化（学界称之为气温直减率），具体的应用需求是建立温度与高程、经纬度的线性模型，并计算模型的F统计量和对应的p值。
数据存在一个叫myData的data.frame中，数据前几列长下面这样：
ID X Y Ele Temp
1 50353 126.6333 51.73333 173.9 0.055024
2 50468 127.4500 50.25000 166.4 1.619074
3 50564 127.3500 49.43333 234.5 1.468917
4 50674 129.4333 48.56667 404.5 0.954324
5 50774 128.8333 47.70000 264.8 2.621272
6 50778 132.5333 47.66667 57.5 4.410369
需要求Temp与Ele、X、Y的线性关系:
myModel = lm(Temp ~ Ele + X + Y, data = myData)
summary(myModel)
Call:
lm(formula = Temp ~ Ele + X + Y, data = myData)
Residuals:
Min 1Q Median 3Q Max
-1.00352 -0.32245 -0.01634 0.39690 1.27065
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 27.8974176 5.0653852 5.507 3.76e-06
Ele -0.0053912 0.0005171 -10.425 3.97e-12
X 0.1243715 0.0423171 2.939 0.00588
Y -0.8286338 0.0377717 -21.938 < 2e-16
(Intercept) ***
Ele ***
X **
Y ***
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.5099 on 34 degrees of freedom
Multiple R-squared: 0.9383, Adjusted R-squared: 0.9328
F-statistic: 172.3 on 3 and 34 DF, p-value: < 2.2e-16
mySummary = summary(myModel)
mySummary[1]或mySummary$call：模型公式
mySummary[3]或mySummary$residuals：模型预测的残差
mySummary[4]或mySummary$coefficient：模型截距，每个变量系数、计算误差、t-value、pr(>t)，放在一个list中，可按需提取
mySummary[6]或mySummary$sigma：residual standard error
mySummary[7]或mySummary$df：自由度
mySummary[8]或mySummary$r.squared：R squared
mySummary[9]或mySummary$adj.r.squared： adjust R squared
mySummary[10]或mySummary$fstatistic：F统计量&一个分子和一个分母值，同样放在list里，可以按需提取
mySummary[11]或mySummarycov.unscaled：变量的相关系数矩阵
mySummary[11]......：NA
这个时候你会发现最后那个p-value不见了，它是根据mySummary$fstatistic里面的几个值求出来的
f = mySummary$fstatistic
p_value = pf(as.numeric(f[1]), as.numeric(f[2]), as.numeric(f[3]), lower.tail = FALSE)
*****************************************************************************************************************
继续气温直减率的问题
lapse_rates = coef(myModel)[2]
f = summary(myModel)$fstatistic
f_statistic = f[1]
p_value = pf(as.numeric(f[1]), as.numeric(f[2]), as.numeric(f[3]), lower.tail = FALSE)
result = c(lapse_rates,f_statistic, p_value)
names(result) = c("lapse_rates","F-statistic","p-value")
print(result)
lapse_rates F-statistic p-value
-5.391219e-03 1.723034e+02 1.261248e-20
有些地方不知道自己写了些啥，他日等我略有长进再来改（2016-8-31 13:05）
