使用R Studio的Rmd文档格式。
完整作业Rmd文档：
---title: "EXP-Assignment-1"output:html_document: defaultword_document: defaultpdf_document: default---```{r setup, include=FALSE}knitr::opts_chunk$set(echo = TRUE)```## 实验作业1 （10%）加载使用的libraries```{r}# Testoptions(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高血压转化为factoroutcome$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)#行标签名，返回是listclass(outcome)#变量类型，结果为："tbl_df""tbl" "data.frame"str(outcome)#查看数据结构，包括行列数、行标签、各变量类型、数据summary(outcome)#查看数据统计信息#variables <-vroom("variables.csv")variables <-vroom(file.choose())#选择文件加载variables$ID <- as.factor(variables$ID)#将数据中分类变量转化为factorvariables$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)#行标签名，返回是liststr(variables)#查看数据结构，包括行列数、行标签、各变量类型、数据summary(variables)#查看数据统计信息```outcome数据统计量：极值包括极大值Max：样本中，数值最大的数据和极小值Min：样本中，数值最小的数据均值Mean：样本中，所有数值的和除以样本数中位数Median：样本中，最中间的一个数值或两个数的均值ID birthyear sysbp hypertension1000553: 1 1944 : 269 Min. : 90.0 0:36651001822: 1 1947 : 269 1st Qu.:125.0 1:13131005779: 1 1946 : 267 Median :136.01006569: 1 1943 : 221 Mean :138.21006871: 1 1948 : 221 3rd Qu.:150.01007611: 1 1945 : 207 Max. :243.0(Other):4972 (Other):3524 NA's :3 variables数据统计量：极值包括极大值Max：样本中，数值最大的数据和极小值Min：样本中，数值最小的数据均值Mean：样本中，所有数值的和除以样本数中位数Median：样本中，最中间的一个数值或两个数的均值alcohol pa1 sleepheightweight 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:12430:1265 0:1265 1st Qu.: 97.0 1:12481:1173 1:1250 Median :102.0 2:12462:1251 2:1280 Mean :103.5 3:12413: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)aprint("缺失值偏少，可以删除，可以用均值,中位数填补")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分）已知患者纳入研究的时间为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#同理，birthmonthint_birthmonth <- as.vector(data_22$birthmonth)class(int_birthmonth)#类型为"character"int_birthmonth <- as.integer(int_birthmonth)#转化为整型，每个数减去最小值了int_birthmonth#data_22$baseline_age <- - int_birthyear - int_birthmonth %/%10#大于9就多算了一岁，减去1#等价于下面data_22 <- data_22 %>% mutate(baseline_age = - 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_levelstr(data_22)```三、一元线性回归（40分） 1.（10分）请选择合适的方法，分析2组变量之间的关系：身高height和吸烟smoking、身高height和体重weight，并用文字解释分析结果；```{r}data_31 = data_new#使用新的变量，以免改变之前的变量#身高height和吸烟smoking是连续变量 & 分类变量（多分类）#将smkoing的factor0、1、2转化为对应的never、ever、currentdata_31$smoking <- factor(data_31$smoking, levels=c(0,1,2), labels=c("never", "ever", "current"))#对于分类变量查看分布my_smoking <- table(data_31$smoking) #频数 my_smokingprop.table(my_smoking) # 百分比# 连续变量是否为正态分布#hist(rnorm(1000))#rnorm(1000)，生成1000满足正态分布的向量，hist将向量画图显示summary(data_31$height)hist(data_31$height)#又：用偏度skewness和峰度Kurtosis检验正态性# 正态分布偏度等于0skewness(data_31$height) # 偏度系数 agostino.test(data_31$height) # 偏度检验# 正态分布的峰度为3 # 当比正态分布的低时，峰度小于3；高，大于3kurtosis(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_smokingnever ever current 2582 1660476 prop.table(my_smoking) # 百分比neverever current 0.5472658 0.3518440 0.1008902summary(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 testdata: data_31$heightskew = 0.18283, z = 5.09332, p-value = 3.519e-07alternative hypothesis: data have a skewnesskurtosis(data_31$height) # 峰度系数 [1] 2.497015接近3，接近满足正态分布anscombe.test(data_31$height) # 峰度检验Anscombe-Glynn kurtosis testdata: data_31$heightkurt = 2.4970, z = -9.5393, p-value < 2.2e-16alternative hypothesis: kurtosis is not equal to 3# ANOVA：单变量方差分析fit <- aov(log(data_31$height) ~ data_31$smoking) model.tables(fit,"means") # 样本均值；每组均值和频数Tables of meansGrand mean5.124178 data_new$smoking neverever current5.12 5.129 5.128rep 2713.00 1726.000 512.000summary(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检验正态性# 正态分布偏度等于0skewness(data_31$weight) # 偏度系数 agostino.test(data_31$weight) # 偏度检验# 正态分布的峰度为3 # 当比正态分布的低时，峰度小于3；高，大于3kurtosis(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检验正态性# 正态分布偏度等于0skewness(data_31$weight) # 偏度系数 [1] 0.664922接近0agostino.test(data_31$weight) # 偏度检验D'Agostino skewness testdata: data_31$weightskew = 0.66492, z = 17.06403, p-value < 2.2e-16alternative hypothesis: data have a skewness# 正态分布的峰度为3 # 当比正态分布的低时，峰度小于3；高，大于3kurtosis(data_31$weight) # 峰度系数 [1] 3.690707大于3anscombe.test(data_31$weight) # 峰度检验Anscombe-Glynn kurtosis testdata: data_31$weightkurt = 3.6907, z = 7.3633, p-value = 1.795e-13alternative 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 correlationdata: data_new$height and data_new$weightt = 45.578, df = 4976, p-value < 2.2e-16alternative hypothesis: true correlation is not equal to 095 percent confidence interval:0.5228007 0.566sample 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 regressionfitsummary(fit) # Check model results###与fit中的lm不同，abline中的lm中，x，y要互换abline(lm(data_32$sysbp~data_32$WaistHipRatio), col = 'blue') # Add fitting lineres <- 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 regressionCall:lm(formula = data_32$WaistHipRatio ~ data_32$sysbp)Coefficients:#回归系数(Intercept) data_32$sysbp 0.727558 0.001029 回归方程：WaistHipRatio = 0.727557819 + 0.001028916 * sysbpconfint(res)#置信区间2.5 %97.5 %(Intercept) 0.7089371098 0.746178529sysbp 0.0008952223 0.00116260997.5%的置信水平在[0.0008952223 0.001162609]anova(res)#残差分析Analysis of Variance TableResponse: WaistHipRatioDf 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 ‘ ’ 1vcov(res)#返回方差和协方差矩阵(Intercept) sysbp(Intercept) 9.021387e-05 -6.419026e-07sysbp -6.419026e-07 4.650518e-09AIC(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 regressionfit1summary(fit1) # Check model results###与fit中的lm不同，abline中的lm中，x，y要互换abline(lm(data_34$hip~data_34$waist), col = 'blue') # Add fitting lineres1 <- 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.406991data_34$hip 1.049244 1.10719197.5%的置信水平在[1.049244 1.107191]#fitted(res1)#拟合#residuals(res1)#残差anova(res1)#残差分析Analysis of Variance TableResponse: data_34$waistDf Sum Sq Mean Sq F value Pr(>F) data_34$hip 1 440466 440466 5322.5 < 2.2e-16 ***Residuals 4716 39027283 ---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1vcov(res1)#返回方差和协方差矩阵(Intercept) data_34$hip(Intercept) 2.35476527 -0.0225942381data_34$hip -0.02259424 0.0002184212AIC(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 regressionfit2summary(fit2) # Check model results#lm中，x，y要互换abline(lm(data_new$BMI~data_new$WaistHipRatio), col = 'blue') # Add fitting lineres2 <- 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.6345840.008574 summary(fit2) # Check model resultsCall: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 ‘ ’ 1Residual standard error: 0.07991 on 4716 degrees of freedomMultiple 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 lineres2 <- with(data_new, fit2)#第一个参数是数据集res2Call:lm(formula = data_new$WaistHipRatio ~ data_new$BMI)Coefficients:(Intercept) data_new$BMI 0.6345840.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 ‘ ’ 1Residual standard error: 0.07991 on 4716 degrees of freedomMultiple R-squared: 0.1996,Adjusted R-squared: 0.1994 F-statistic: 1176 on 1 and 4716 DF, p-value: < 2.2e-16coefficients(res2)#回归系数(Intercept) data_new$BMI 0.634584117 0.008573744 confint(res2)#置信区间2.5 %97.5 %(Intercept) 0.620956367 0.64821187data_new$BMI 0.008083548 0.00906394#fitted(res2)#拟合#res2iduals(res2)#残差anova(res2)#残差分析Analysis of Variance TableResponse: data_new$WaistHipRatioDf 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 ‘ ’ 1vcov(res2)#返回方差和协方差矩阵(Intercept) data_new$BMI(Intercept) 4.832025e-05 -1.713582e-06data_new$BMI -1.713582e-06 6.25e-08AIC(res2) #AIC即赤池信息准则，用于判断模型的拟合优度，越小越好-10449.77回归方程：WaistHipRatio = 0.634584117 + 0.008573744 * BMI两者是线性关系，绝大多数点在回归曲线附近有残差图可知，曲线有些弯曲，可能未对异常值进行处理，导致模型拟合的不是很好。选择的模型不是很好，需要进行修正。
