R语言数据读取、清洗、一元线性回归

稳健的不高冷的强哥
已于 2022-04-21 12:15:22 修改
阅读量1k
点赞数
分类专栏： R语言文章标签： r语言机器学习
于 2022-04-21 12:09:46 首次发布
本文链接：https://blog.csdn.net/qq_40925617/article/details/124318395
版权
R语言专栏收录该内容
1 篇文章 0 订阅
订阅专栏
使用R Studio的Rmd文档格式。
完整作业Rmd文档：
---
title: "EXP-Assignment-1"
output:
  html_document: default
  word_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## 实验作业1 （10%）

加载使用的libraries
```{r}
# Test
options(download.file.method="libcurl")
p <- c("dplyr", "vroom", "vioplot", "tidyverse", "moments", "shiny")#批量下载并导入包
#installed.packages()[, "Package"] #查看已经下载的包

ipak <- function(x){
    new.x <- x[!(x %in% installed.packages()[, "Package"])]
    #print(new.x)#打印输出一下，未下载的包
    if (length(new.x))
        install.packages(new.x, dependencies = TRUE)
    sapply(x, require, character.only = TRUE)
}
ipak(p)

```
(test: 文字解释)

一、数据框的基本操作（30分）  
1.（10分）请分别对两个数据集进行读取和数据概览（行列数、变量名、变量类型），并用*文字*进行简单统计描述（极值、均值、中位数）；
```{r}
#outcome <- read.csv("outcome.csv")
# Recommend vroom, quick speed for large data file, 
# automatically identified delimiters, column names and data types
#outcome <- vroom("outcome.csv")
outcome <- vroom(file.choose())#选择文件加载
outcome$ID <- as.factor(outcome$ID)#将ID、birthday高血压转化为factor
outcome$birthyear <- as.factor(outcome$birthyear)
outcome$hypertension <- as.factor(outcome$hypertension)

outcome <- distinct(outcome, ID, .keep_all = T)#根据ID，去重,其他保留
dim(outcome)#行列数
colnames(outcome)#行标签名
names(outcome)#行标签名，返回是list
class(outcome)#变量类型，结果为："tbl_df"     "tbl"        "data.frame"
str(outcome)#查看数据结构，包括行列数、行标签、各变量类型、数据
summary(outcome)#查看数据统计信息

#variables <-vroom("variables.csv")
variables <-vroom(file.choose())#选择文件加载
variables$ID <- as.factor(variables$ID)#将数据中分类变量转化为factor
variables$birthyear <- as.factor(variables$birthyear)
variables$birthmonth <- as.factor(variables$birthmonth)
variables$gender <- as.factor(variables$gender)
variables$race <- as.factor(variables$race)
variables$education <- as.factor(variables$education)
variables$smoking <- as.factor(variables$smoking)
variables$alcohol <- as.factor(variables$alcohol)
variables$pa1 <- as.factor(variables$pa1)
variables$sleep <- as.factor(variables$sleep)
variables$weightrank <- as.factor(variables$weightrank)
variables$waistrank <- as.factor(variables$waistrank)
variables$hiprank <- as.factor(variables$hiprank)

variables <- distinct(variables, ID, .keep_all = T)#根据ID，去重
dim(variables)#行列数
colnames(variables)#行标签名
class(variables)#变量类型
names(variables)#行标签名，返回是list
str(variables)#查看数据结构，包括行列数、行标签、各变量类型、数据
summary(variables)#查看数据统计信息
```
 outcome数据统计量：
 极值包括极大值Max：样本中，数值最大的数据
 和极小值Min：样本中，数值最小的数据
 均值Mean：样本中，所有数值的和除以样本数
 中位数Median：样本中，最中间的一个数值或两个数的均值
        ID         birthyear        sysbp       hypertension
 1000553:   1   1944   : 269   Min.   : 90.0   0:3665      
 1001822:   1   1947   : 269   1st Qu.:125.0   1:1313      
 1005779:   1   1946   : 267   Median :136.0               
 1006569:   1   1943   : 221   Mean   :138.2               
 1006871:   1   1948   : 221   3rd Qu.:150.0               
 1007611:   1   1945   : 207   Max.   :243.0               
 (Other):4972   (Other):3524   NA's   :3 
 
 variables数据统计量：
 极值包括极大值Max：样本中，数值最大的数据
 和极小值Min：样本中，数值最小的数据
 均值Mean：样本中，所有数值的和除以样本数
 中位数Median：样本中，最中间的一个数值或两个数的均值
 alcohol       pa1        sleep          height          weight           waist       
 0   : 239   0   :2098   0   :3376   Min.   :142.0   Min.   : 38.60   Min.   : 57.00  
 1   : 203   1   :1184   1   :1178   1st Qu.:161.0   1st Qu.: 66.40   1st Qu.: 80.00  
 2   :4515   2   :1586   2   : 372   Median :168.0   Median : 76.60   Median : 90.00  
 NA's:  21   NA's: 110   NA's:  52   Mean   :168.3   Mean   : 77.99   Mean   : 90.27  
                                     3rd Qu.:175.0   3rd Qu.: 87.57   3rd Qu.: 99.00  
                                     Max.   :198.0   Max.   :172.20   Max.   :154.00  
                                                                                      
      hip        weightrank waistrank hiprank 
 Min.   : 78.0   0:1243     0:1265    0:1265  
 1st Qu.: 97.0   1:1248     1:1173    1:1250  
 Median :102.0   2:1246     2:1251    2:1280  
 Mean   :103.5   3:1241     3:1289    3:1183  
 3rd Qu.:108.0                                
 Max.   :154.0 
2.（10分）请将两个数据集分别按行、列进行直接组合，如果遇到无法组合的Error报错，请解释报错原因；
```{r}
#按行合并
#rdata <- rbind(outcome, variables)#行数不对，报错

#按列合并
cdata <- cbind(outcome, variables)#ID和birthyear重复，列合并不好
#推荐下面几种方法合并
#data2 <- merge(outcome, variables, by="ID")
#data2 <- full_join(outcome,variables)#列合并，自动将相同的列合并
#left_join
#inner_join
```

3.（10分）请将两个数据集按ID合并为新的数据集，命名为data2，并对data2进行第1题中的数据概览和描述。
```{r}
data2 <- merge(outcome, variables, by="ID")#按ID合并，但会出现两个birthyear
#data2 <- full_join(outcome, variables, by='ID')
dim(data2)#行列数
colnames(data2)#行标签名
names(data2)#行标签名
class(data2)#变量类型
str(data2)#变量内部结构，行标签，变量类型，数据
summary(data2)#总结统计量
```
二、常用的数据处理（30分）  
1.（10分）数据集data2是否存在缺失值？如果存在，请选择你认为合适的方式处理缺失值，并将处理后的数据集命名为data_new；
```{r}
data2 <- merge(outcome,variables)#按相同的数据标签自动合并
#查看变量缺失值情况，无缺失值为100 #apply函数，1为行，对行操作，2为列，对列操作
apply(data2, 2, function(x){sum(is.na(x))/length(x)*100})
print("查看行缺失情况，1为无确实行")
#complete.cases(data2)，查看行是否完整
a <- sum(!complete.cases(data2)) / nrow(data2)
a
print("缺失值偏少，可以删除，可以用均值,中位数填补")
if(a <= 0.1){
  data_new <- na.omit(data2)#缺失行比例小于0.05,直接删除
}else{#使用均值填补，
  data_new <- data2
  #此数据种类较多，均值填补比较麻烦，
  #apply(data2,2,function(x){x[is.na(x)] = mean(x, na.rm = TRUE)})
}
summary(data_new)
```
2.（10分）已知患者纳入研究的时间为2008年9月，请计算每个患者的基线年龄，并在data_new中添加一列，命名为baseline_age；
```{r}
data_22 = data_new#使用新的变量，以免改变之前的变量
#将birthyear和birthmonth转化为整数
#factor直接转化为整型，结果不对，貌似每个数减去了某个数
#因此先factor转化为向量"character"，再转化为整型
int_birthyear <- as.vector(data_22$birthyear)
class(int_birthyear)#类型为"character"
int_birthyear <- as.integer(int_birthyear)#转化为整型，每个数减去最小值了
int_birthyear
#同理，birthmonth
int_birthmonth <- as.vector(data_22$birthmonth)
class(int_birthmonth)#类型为"character"
int_birthmonth <- as.integer(int_birthmonth)#转化为整型，每个数减去最小值了
int_birthmonth

#data_22$baseline_age <- 2008 - int_birthyear - int_birthmonth %/%10#大于9就多算了一岁，减去1
#等价于下面
data_22 <- data_22 %>% mutate(baseline_age = 2008 - int_birthyear - int_birthmonth%/%9, .keep="all") 
str(data_22)#查看信息
```
3.（10分）请将基线年龄变量baseline_age分层赋值，并在data_new中添加一列，命名为age_level。（分层赋值要求：0-20岁：1， 21-40岁：2， 41-50岁：3， 51-60岁：4， 61-70岁：5， 70岁以上：6）
```{r}
head(data_22$baseline_age)#查看以下数据
age_level <- dplyr::case_when(#此处&与&&结果不同,必须使用&
  data_22$baseline_age <= 20 ~ 1, #~ 赋值的意思
  (data_22$baseline_age > 20 & data_22$baseline_age <= 40) ~ 2,
  (data_22$baseline_age > 40 & data_22$baseline_age <= 50) ~ 3,
  (data_22$baseline_age > 50 & data_22$baseline_age <= 60) ~ 4,
  (data_22$baseline_age > 60 & data_22$baseline_age <= 70) ~ 5,
  (data_22$baseline_age >70) ~ 6,
)
#str(age_level)
data_22$age_level = age_level
str(data_22)
```
三、一元线性回归（40分）  
1.（10分）请选择合适的方法，分析2组变量之间的关系：身高height和吸烟smoking、身高height和体重weight，并用文字解释分析结果；
```{r}
data_31 = data_new#使用新的变量，以免改变之前的变量
#身高height和吸烟smoking是连续变量 & 分类变量（多分类）
#将smkoing的factor0、1、2转化为对应的never、ever、current
data_31$smoking <- factor(data_31$smoking, levels=c(0,1,2), labels=c("never", "ever", "current"))

#对于分类变量查看分布
my_smoking <- table(data_31$smoking) #频数 
my_smoking
prop.table(my_smoking) # 百分比

# 连续变量是否为正态分布
#hist(rnorm(1000))#rnorm(1000)，生成1000满足正态分布的向量，hist将向量画图显示
summary(data_31$height)
hist(data_31$height)
#又：用偏度skewness和峰度Kurtosis检验正态性
# 正态分布偏度等于0
skewness(data_31$height) # 偏度系数 
agostino.test(data_31$height)  # 偏度检验

# 正态分布的峰度为3 
# 当比正态分布的低时，峰度小于3；高，大于3
kurtosis(data_31$height) # 峰度系数 
anscombe.test(data_31$height)  # 峰度检验

# ANOVA：单变量方差分析
fit <- aov(log(data_31$height) ~ data_31$smoking) 
model.tables(fit,"means") # 样本均值；每组均值和频数
summary(fit) #F统计量和p值
```
height是连续性变量，smoking是多分类型变量
my_smoking
  never    ever current 
   2582    1660     476 
prop.table(my_smoking) # 百分比
    never      ever   current 
0.5472658 0.3518440 0.1008902

summary(data_31$height)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  142.0   161.5   168.0   168.4   175.0   198.0 
  
skewness(data_31$height) # 偏度系数 
[1] 0.1828344接近0，接近满足正态分布

agostino.test(data_31$height)  # 偏度检验
	D'Agostino skewness test
data:  data_31$height
skew = 0.18283, z = 5.09332, p-value = 3.519e-07
alternative hypothesis: data have a skewness

kurtosis(data_31$height) # 峰度系数 
[1] 2.497015接近3，接近满足正态分布

anscombe.test(data_31$height)  # 峰度检验
	Anscombe-Glynn kurtosis test
data:  data_31$height
kurt = 2.4970, z = -9.5393, p-value < 2.2e-16
alternative hypothesis: kurtosis is not equal to 3

# ANOVA：单变量方差分析
fit <- aov(log(data_31$height) ~ data_31$smoking) 
model.tables(fit,"means") # 样本均值；每组均值和频数
Tables of means
Grand mean
5.124178 
 data_new$smoking 
      never     ever current
       5.12    5.129   5.128
rep 2713.00 1726.000 512.000

summary(fit) #F统计量和p值
                   Df Sum Sq Mean Sq F value   Pr(>F)    
data_new$smoking    2  0.095 0.04767   16.22 9.55e-08 ***
Residuals        4948 14.545 0.00294
```{r}
# 身高height和体重weight是连续变量 & 连续变量
#使用相关分析correlation
# 连续变量是否为正态分布
summary(data_new$height)#身高
hist(data_new$height)
summary(data_new$weight)#体重
hist(data_new$weight)

#又：用偏度skewness和峰度Kurtosis检验正态性
# 正态分布偏度等于0
skewness(data_31$weight) # 偏度系数 
agostino.test(data_31$weight)  # 偏度检验

# 正态分布的峰度为3 
# 当比正态分布的低时，峰度小于3；高，大于3
kurtosis(data_31$weight) # 峰度系数 
anscombe.test(data_31$weight)  # 峰度检验

# 相关系数的计算,两变量均服从正态分布 pearson相关系数 
cor(data_new$height, data_new$weight, method="pearson") #求相关系数
# 求相关系数的P值
cor.test(data_new$height, data_new$weight, method="pearson") 

```
height和weight都是连续性变量
 #使用相关分析correlation
 # 连续变量是否为正态分布
 summary(data_new$height)#身高
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  142.0   161.5   168.0   168.4   175.0   198.0 
 hist(data_new$height)
 summary(data_new$weight)#体重
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  39.90   66.40   76.60   77.95   87.50  172.20 
 hist(data_new$weight)
 
 #又：用偏度skewness和峰度Kurtosis检验正态性
 # 正态分布偏度等于0
 skewness(data_31$weight) # 偏度系数 
[1] 0.664922接近0

 agostino.test(data_31$weight)  # 偏度检验
	D'Agostino skewness test
data:  data_31$weight
skew = 0.66492, z = 17.06403, p-value < 2.2e-16
alternative hypothesis: data have a skewness

 # 正态分布的峰度为3 
 # 当比正态分布的低时，峰度小于3；高，大于3
 kurtosis(data_31$weight) # 峰度系数 
[1] 3.690707大于3

 anscombe.test(data_31$weight)  # 峰度检验
	Anscombe-Glynn kurtosis test
data:  data_31$weight
kurt = 3.6907, z = 7.3633, p-value = 1.795e-13
alternative hypothesis: kurtosis is not equal to 3

 相关系数的计算,两变量均服从正态分布 pearson相关系数 
 cor(data_new$height, data_new$weight, method="pearson") #求相关系数
[1] 0.5404557
求相关系数的P值
 cor.test(data_new$height, data_new$weight, method="pearson") 
	Pearson's product-moment correlation
data:  data_new$height and data_new$weight
t = 45.578, df = 4976, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5228007 0.5620066
sample estimates:
      cor 
0.5426992 
相关系数：r=0.5426992 > 0, 接近1，说明两者相关
p值：p<2.2e-16<0.001, 相关系数有很强的确定性

2.（10分）请拟合一元线性回归模型，分析腰臀比(Waist-to-Hip Ratio, WHR=waist/hip)和收缩压(sysbp)的关系，并绘制散点图和回归线，用*文字*解释拟合出的模型结果；
```{r}
data_32 = data_new#使用新的变量，以免改变之前的变量
data_32$WaistHipRatio <- data_32$waist / data_32$hip#计算腰臀比
#通过散点图初步判断变量之间的关联
plot(data_32$WaistHipRatio, data_32$sysbp, type="p", xlab="WaistHipRatio", ylab="sysbp")
#初步进行线性回归
# with()函数指定数据集
fit <- lm(data_32$WaistHipRatio~data_32$sysbp) # Linear regression
fit
summary(fit) # Check model results
###与fit中的lm不同，abline中的lm中，x，y要互换
abline(lm(data_32$sysbp~data_32$WaistHipRatio), col = 'blue') # Add fitting line

res <- with(data_32, fit)#第一个参数是数据集
res

#模型诊断
par(mfrow=c(2,2))#将图画在2x2矩阵上
plot(res)

#模型解读 
summary(res)
coefficients(res)#回归系数
confint(res)#置信区间
#fitted(res)#拟合
#residuals(res)#残差
anova(res)#残差分析
vcov(res)#返回方差和协方差矩阵
AIC(res)  #AIC即赤池信息准则，用于判断模型的拟合优度，越小越好

#对原数据每个进行预测
result <- predict(res)
#result#数据太多，因此不输出
```
散点图：无特殊结构；表示符合线性
残差图：残差与真实值之间的关系画图，判断变量是否呈线性，线性，则残差与回归预测值不相关，
Q-Q图：判断残差正态分布，图中散点都集中在QQ图直线上，说明残差满足正态分布，验证数据服从正态分布
Scale-Location图：检查等方差假设，判断残差方差齐性：散点图显示两个变量无关联，表示残差方差齐性
residual vs leverage: 用于检查分析数据中是否存在极端值？？？Cook距离。咱们在线性模型里用Cook距离分析一个点是否很是“influential。”通常来讲距离大于0.5的点就须要引发注意了，需注意，即便R将这些特殊的点标记了出来，也不等于他们必定须要被删除。仍是要参考Cook距离的绝对大小
fit <- lm(data_32$WaistHipRatio~data_32$sysbp) # Linear regression
Call:
lm(formula = data_32$WaistHipRatio ~ data_32$sysbp)

Coefficients:#回归系数
  (Intercept)  data_32$sysbp  
     0.727558       0.001029  
回归方程：WaistHipRatio = 0.727557819 + 0.001028916 * sysbp

confint(res)#置信区间
                   2.5 %      97.5 %
(Intercept) 0.7089371098 0.746178529
sysbp       0.0008952223 0.001162609
97.5%的置信水平在[0.0008952223 0.001162609]

anova(res)#残差分析
Analysis of Variance Table
Response: WaistHipRatio
            Df Sum Sq Mean Sq F value    Pr(>F)    
sysbp        1  1.733 1.73263  227.65 < 2.2e-16 ***
Residuals 4716 35.894 0.00761      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

vcov(res)#返回方差和协方差矩阵
              (Intercept)         sysbp
(Intercept)  9.021387e-05 -6.419026e-07
sysbp       -6.419026e-07  4.650518e-09

AIC(res)  #AIC即赤池信息准则，用于判断模型的拟合优度，越小越好
-9621.985
回归方程：
WaistHipRatio = 0.727557819 + 0.001028916 * sysbp
两者是线性关系，绝大多数点在回归曲线附近
3.（8分）请绘制第2题中模型残差的散点图，并分析模型可能存在的问题；
```{r}
#残差散点图
plot(res,which = 1)
```
残差图：残差与真实值之间的关系画图，判断变量是否呈线性，线性，则残差与回归预测值不相关，
曲线有些弯曲，不是线性关系，说明残差和估计值有关，模型效果不是很好。
未对异常值进行处理，导致模型拟合的不是很好。
选择的模型不是很好，需要进行修正。

4.（12分）请任选2个你感兴趣的变量作为自变量和因变量，描述你的研究问题，并进行线性回归分析。（研究问题的选择-2分、模型拟合和诊断-4分、结果解释和可视化正确且充实-6分）
```{r}
data_34 = data_new#使用新的变量，以免改变之前的变量
#线性回归：针对自变量和因变量都是连续性变量的分析 
#数据集中，只有sysbp、height、weight、waist、hip几个连续性变量
#研究身腰Waist和臀waist的关系
#类似第三题的分析
#通过散点图初步判断变量之间的关联，由散点图可知，Waist和臀waist呈现线性关系
plot(data_34$waist, data_34$hip, type="p", xlab="waist", ylab="hip")
#初步进行线性回归
# with()函数指定数据集
fit1 <- lm(data_34$waist~data_34$hip) # Linear regression
fit1
summary(fit1) # Check model results
###与fit中的lm不同，abline中的lm中，x，y要互换
abline(lm(data_34$hip~data_34$waist), col = 'blue') # Add fitting line

res1 <- with(data_34, fit1)#第一个参数是数据集
res1

#模型诊断
par(mfrow=c(2,2))#将图画在2x2矩阵上
plot(res1)

#模型解读 
summary(res1)
coefficients(res1)#回归系数
confint(res1)#置信区间
#fitted(res1)#拟合
#residuals(res1)#残差
anova(res1)#残差分析
vcov(res1)#返回方差和协方差矩阵
AIC(res1)  #AIC即赤池信息准则，用于判断模型的拟合优度，越小越好

#对原数据每个进行预测
result1 <- predict(res1)
#result1#数据太多，因此不输出
```
散点图：无特殊结构；表示符合线性，由散点图可知，Waist和臀waist呈现线性关系
残差图：残差与真实值之间的关系画图，判断变量是否呈线性，线性，则残差与回归预测值不相关，
Q-Q图：判断残差正态分布，图中散点都集中在QQ图直线上，说明残差满足正态分布，验证数据服从正态分布
Scale-Location图：检查等方差假设，判断残差方差齐性：散点图显示两个变量无关联，表示残差方差齐性
residual vs leverage: 用于检查分析数据中是否存在极端值？？？Cook距离。咱们在线性模型里用Cook距离分析一个点是否很是“influential。”通常来讲距离大于0.5的点就须要引发注意了，需注意，即便R将这些特殊的点标记了出来，也不等于他们必定须要被删除。仍是要参考Cook距离的绝对大小

 coefficients(res1)#回归系数
(Intercept) data_34$hip 
 -21.415376    1.078218 
 confint(res1)#置信区间
                 2.5 %     97.5 %
(Intercept) -24.423761 -18.406991
data_34$hip   1.049244   1.107191
97.5%的置信水平在[1.049244   1.107191]

 #fitted(res1)#拟合
 #residuals(res1)#残差
 anova(res1)#残差分析
Analysis of Variance Table

Response: data_34$waist
              Df Sum Sq Mean Sq F value    Pr(>F)    
data_34$hip    1 440466  440466  5322.5 < 2.2e-16 ***
Residuals   4716 390272      83                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 vcov(res1)#返回方差和协方差矩阵
            (Intercept)   data_34$hip
(Intercept)  2.35476527 -0.0225942381
data_34$hip -0.02259424  0.0002184212
 AIC(res1)  #AIC即赤池信息准则，用于判断模型的拟合优度，越小越好
[1] 34227.25

回归方程：
waist = -21.415376 + 1.078218 * hip
两者是线性关系，绝大多数点在回归曲线附近
由残差图可知，曲线有些弯曲，可能未对异常值进行处理，导致模型拟合的不是很好。
选择的模型不是很好，需要进行修正。
```{r}

#研究身高体重BMI和腰臀比(Waist-to-Hip Ratio, WHR=waist/hip)的关系
#BMI <- data_new %>% mutate(BMI = round(weight/(height*1.0/100)^2))
#BMI %>% ggplot() + aes(x = age, y = BMI) + geom_point() + geom_smooth(method='lm')
#类似第三题的分析
data_new$BMI <- round(data_new$weight / (data_new$height*1.0/100)^2, digits = 1)#四舍五入，一位小数
data_new$WaistHipRatio <- data_new$waist / data_new$hip#计算腰臀比
#summary(data_new)#查看添加信息
#通过散点图初步判断变量之间的关联
plot(data_new$WaistHipRatio, data_new$BMI, type="p", xlab="WaistHipRatio", ylab="BMI")
#初步进行线性回归
# with()函数指定数据集
fit2 <- lm(data_new$WaistHipRatio~data_new$BMI)#Linear regression
fit2
summary(fit2) # Check model results
#lm中，x，y要互换
abline(lm(data_new$BMI~data_new$WaistHipRatio), col = 'blue') # Add fitting line

res2 <- with(data_new, fit2)#第一个参数是数据集
res2

#模型诊断
par(mfrow=c(2,2))#将图画在2x2矩阵上
plot(res2)

#模型解读 
summary(res2)
coefficients(res2)#回归系数
confint(res2)#置信区间
#fitted(res2)#拟合
#res2iduals(res2)#残差
anova(res2)#残差分析
vcov(res2)#返回方差和协方差矩阵
AIC(res2)  #AIC即赤池信息准则，用于判断模型的拟合优度，越小越好

#对原数据每个进行预测
result <- predict(res2)
#result#数据太多，因此不输出
```


散点图：无特殊结构；表示符合线性
残差图：残差与真实值之间的关系画图，判断变量是否呈线性，线性，则残差与回归预测值不相关，
Q-Q图：判断残差正态分布，图中散点都集中在QQ图直线上，说明残差满足正态分布，验证数据服从正态分布
Scale-Location图：检查等方差假设，判断残差方差齐性：散点图显示两个变量无关联，表示残差方差齐性
residual vs leverage: 用于检查分析数据中是否存在极端值？？？Cook距离。咱们在线性模型里用Cook距离分析一个点是否很是“influential。”通常来讲距离大于0.5的点就须要引发注意了，需注意，即便R将这些特殊的点标记了出来，也不等于他们必定须要被删除。仍是要参考Cook距离的绝对大小

Call:
lm(formula = data_new$WaistHipRatio ~ data_new$BMI)

Coefficients:
 (Intercept)  data_new$BMI  
    0.634584      0.008574  

 summary(fit2) # Check model results

Call:
lm(formula = data_new$WaistHipRatio ~ data_new$BMI)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.293102 -0.059664  0.003307  0.060241  0.239345 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.634584   0.006951   91.29   <2e-16 ***
data_new$BMI 0.008574   0.000250   34.29   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.07991 on 4716 degrees of freedom
Multiple R-squared:  0.1996,	Adjusted R-squared:  0.1994 
F-statistic:  1176 on 1 and 4716 DF,  p-value: < 2.2e-16

 #lm中，x，y要互换
 abline(lm(data_new$BMI~data_new$WaistHipRatio), col = 'blue') # Add fitting line
 
 res2 <- with(data_new, fit2)#第一个参数是数据集
 res2

Call:
lm(formula = data_new$WaistHipRatio ~ data_new$BMI)

Coefficients:
 (Intercept)  data_new$BMI  
    0.634584      0.008574  

 #模型诊断
 par(mfrow=c(2,2))#将图画在2x2矩阵上
 plot(res2)
 
 #模型解读 
 summary(res2)

Call:
lm(formula = data_new$WaistHipRatio ~ data_new$BMI)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.293102 -0.059664  0.003307  0.060241  0.239345 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.634584   0.006951   91.29   <2e-16 ***
data_new$BMI 0.008574   0.000250   34.29   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.07991 on 4716 degrees of freedom
Multiple R-squared:  0.1996,	Adjusted R-squared:  0.1994 
F-statistic:  1176 on 1 and 4716 DF,  p-value: < 2.2e-16

coefficients(res2)#回归系数
 (Intercept) data_new$BMI 
 0.634584117  0.008573744 
confint(res2)#置信区间
                   2.5 %     97.5 %
(Intercept)  0.620956367 0.64821187
data_new$BMI 0.008083548 0.00906394
#fitted(res2)#拟合
#res2iduals(res2)#残差
anova(res2)#残差分析
Analysis of Variance Table

Response: data_new$WaistHipRatio
               Df  Sum Sq Mean Sq F value    Pr(>F)    
data_new$BMI    1  7.5088  7.5088  1175.8 < 2.2e-16 ***
Residuals    4716 30.1179  0.0064                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
vcov(res2)#返回方差和协方差矩阵
               (Intercept)  data_new$BMI
(Intercept)   4.832025e-05 -1.713582e-06
data_new$BMI -1.713582e-06  6.252020e-08
AIC(res2)  #AIC即赤池信息准则，用于判断模型的拟合优度，越小越好
-10449.77
回归方程：
WaistHipRatio = 0.634584117 + 0.008573744  * BMI
两者是线性关系，绝大多数点在回归曲线附近
有残差图可知，曲线有些弯曲，可能未对异常值进行处理，导致模型拟合的不是很好。
选择的模型不是很好，需要进行修正。