###第1章：基础入门
##第1章-3.R扩展包的10个问题
#安装某个扩展包
#也可以在右下角的框中“packages”栏中的“install”项目中操作
install.packages("devtools")
#查看已安装的扩展包
#也可以在右下角的框中“packages”栏中查看
installed.packages()
#更新：第一行查看未更新的扩展包；第二行执行更新指令。
#也可以在右下角的框中“packages”栏中的“update”项目中操作
old.packages()
update.packages()
#删除拓展包
#也可以在右下角的框中“packages”栏的“×”操作
remove.packages("devtools")
#加载扩展包
#也可以在右下角的框中“packages”栏中相应的扩展包前打√
#拓展包只能一个一个加载，无法同时加载多个拓展包
library("devtools")
#取消加载
#也可以在右下角的框中“packages”栏中相应的扩展包前去掉√
detach("package:devtools", unload=TRUE)
#查找某个扩展包的作用
#在右下角的框中“help”栏中搜索相应扩展包名字
??devtools
#查找某个扩展包的某功能的介绍
#在右下角的框中“help”栏中搜索某个扩展包的功能名
#下面以查找dplyr扩展包的mutate功能演示
install.packages("dplyr")
library("dplyr")
help(mutate,package = "dplyr")
#如何找到我们要用哪种扩展包分析？
#方法一：找相似的文献
#方法二：CRAN Task Views网站：https://cran.r-project.org/web/views/
#方法三：Bioconductor网站：https://www.bioconductor.org/
###第二章：数据的查看与整理
##第二章-1.数据的基本类型和保存形式
#数字的基本类型：
#数字（double/numeric，缩写num）
#整数（integer，缩写int)：1L，2L…
#虚数（complex，缩写cplx）：1+2i，3-5i…
#逻辑（logic，缩写logi）：True，False
#文字（character，缩写chr）
#数字运算需要先安装并加载base扩展包
install.packages("base")
library("base")
#向量(vector): 维度1，数据类型一致
#生成一行包含“1,2.5,4.5”的数字向量
## “<-”表示赋值
dbl_var<-c(1,2.5,4.5)
dbl_var
#生成一行包含“Hello、world”的文字向量
chr_var<-c("Hello",'world')
chr_var
#生成逻辑变量，TRUE和T都可以表示对，FALSE和F都可以表示错
log_var<-c(TRUE,FALSE,T,F)
log_var
#数据类型不一致时，R会自动把向量保存成较为灵活的形式
#eg1：输入实数和整数，输出时都保存成实数形式
dbl_var1<-c(1.5,1L)
dbl_var1
#eg2：输入整数和文本，输出是都保存成文本形式
dbl_var2<-c(2,'a')
dbl_var2
#列表(list):维度1，数据类型可以不一致
x<-list(1:3,"a",c(TRUE, FALSE, TRUE),list(2.3,5.9))
x
#如何快速查看数据结构? 这个功能对于数据量大、结构复杂的列表非常有帮助
str(x)
#矩阵(matrix):维度2，数据类型一致
#赋值：构建了一个宽为3列，长为2行的矩阵，并在其中填充数字1-6
mat_a <- matrix(1:6,ncol=3,nrow=2)
mat_a
#换一种说法，矩阵可以看做是向量构成的
#创建一个向量b，包含6个数字：1-6
mat_b <- 1:6
#给向量添加一个3行2列的矩阵框架
dim(mat_b)<- c(3,2)
#查看
mat_b
#数据框(data.frame):维度2，数据类型不一定一致
#数据框不用设置长和宽
#但是当数值长度不一样时，会出现下例中的z自动填充一样的数值0，但更多情况是直接报错
install.packages("graphics")
library("graphics")
#赋值
df <- data.frame(x=1:3,y=c("a","b","c"),z=0)
#查看
#数据框也可以通过点击右上角框内“Environment”栏目下方的蓝色三角查看代码计算结果
#点解缩小的网格查看表格
str(df)
#阵列结构一般用不到，需要请另行自学
#数据框转变矩阵再转变会数据框
#矩阵转换为向量：as.vector()
#矩阵转换为数据框：as.data.frame()
#数据框转换为矩阵：as.matrix()
#生成数据框并查看
df1<-data.frame(x=1:3,y=c("a","b","c"),z=0)
df1
str(df1)
#转换成矩阵并查看
#发现生成了3x3的矩阵，但是所有格式都变成了文本形式，因为矩阵要求格式一致
mat_df1<-as.matrix(df1)
str(mat_df1)
#再转换回数据框
#重新变回3x3的数据框，但是所有格式都变成了文本形式
new_df1<-as.data.frame(mat_df1)
str(new_df1)
#分类变量：Factor是分类变量的常见形式
#R的数据框中常常把文字保存成分类变量，包含两个部分：数字（value）和文字（level）
#输入性别数据，生成文本形式的分类变量
sex<-c('male','female','male','female','male','male')
str(sex)
#把文本转换成数字形式的分类变量（factor）
sex_factor<-factor(sex,levels=c('male','female'))
str(sex_factor)
##第二章-2.数据的基本运算和数据提取
##第二章-2.数据的基本运算
#数字运算：
#加
1+1
#减
2-1
#乘
2*3
#除
6/2
#次方
2^3
#开方，以下为几的2次方为9
sqrt(9)
#指数函数：表示以e为底数的指数，以下表示e的2次方
exp(2)
#以e为底的自然对数，以下表示求e的几次方为7.390
log(7.390)
#向量运算：
#产生包含1,2,3三个数的向量，查看
vec_a<-c(1,2,3)
vec_a
#向量中各个数字加一
vec_a+1
#向量中各个数字与自己相加
vec_a+vec_a
#向量中各个数字乘以3
vec_a*3
#求向量中各个数字的指数函数
exp(vec_a)
#矩阵运算：
#生成2行3列的矩阵，软件会按顺序把数字先填第一列，再填第二列，再...
mat_a<-matrix(data=c(1,4,4,9,16,25),ncol=3,nrow=2)
mat_a
#矩阵的各个数字加5
mat_a+5
#矩阵的各个数字开方
sqrt(mat_a)
#把矩阵的行变成列，把列变成成行
t(mat_a)
#计算矩阵的内积（没搞懂，要用再说）
mat_a%*%t(mat_a)
##第二章-2.数据的提取
#数据的提取：
#①向量数据的提取
vec_b<-c("a","b","c")
vec_b
#提取向量第1个数据
vec_b[1]
#提取向量第2个数据
vec_b[2]
#提取向量第3个数据
vec_b[3]
#提取向量第1、第3个数据
vec_b[c(1,3)]
#②矩阵数据的提取
mat_b<-matrix(data=c("a","b","c","d"),ncol=2,nrow=2)
mat_b
#提取第1行，第2列的数据
mat_b[1,2]
#提取第2行，第2列的数据
mat_b[2,2]
#提取第1行所有的数据
mat_b[1,]
#提取第2列所有的数据
mat_b[,2]
#③数据框数据的提取
#x、y、z分别为第1、2、3列的数据框；
#x列的数据为1、2、3；y列的数据为a、b、c；z列数据为0、0、0；
#提取方法1:数据框名【所在行数，所在列数】
#提取方法2:数据框名【所在行数，所在列的名字】
df2<-data.frame(x=1:3,y=c("a","b","c"),z=0)
df2
str(df2)
df2[1,1]#显示第一行，第一列（x列）
df2[1,'x']#显示第一行，第一列（x列）
df2[1,]#显示第一行
df2[,'y']#显示第2列（y列）
df2[,1]#显示第1列（x列）
df2$y#取出y中所有数据
##第二章-3.数据的读入与保存
#路径的识别与更改：
#①查看工作路径（资料存储、运行的位置）
#建议每开一个新项目，就创建一个文件夹，把该项目的工作路径都改到该文件夹：
getwd()
#②更改工作路径：
#方法1：代码
#方法2：右下角的框→“file"→“More"→“go to working derectory”→打开相应的文件夹→“More"→“set as Working Derectory”
#方法3：右下角的框→“…”→打开某文件夹→设定该路径为工作路径
setwd("路径名")
#③查看路径是否已经更改：右下角的框→“file"→“More"→“go to working derectory”→房子后面那一堆就是现在的工作路径
#Excel文件的读取
#方法一：代码
#Excel文件在医咖会课程下载
#课程有3种拓展包可以用于读取Excel文件：gdata包、XLConnect包、xslx包（都不行）
#系统自带的readxl包、上网找的rio包：均可行
#自己上网找了一个可行的：rio包
#操作前需要先把工作路径改到有操作Excel文件的文件夹
#readxl包：
install.packages("readxl")
library("readxl")
example1 <- read_excel("YKH.example1.xlsx")
View(example1)
#rio包：
install.packages("rio")
library("rio")
example2<-rio::import("YKH.example1.xlsx")
View(example2)
#Excel文件的读取：
#方法二：鼠标操作
#最上方的"file"→"import dataset"→"from excel"→"browse"→选择Excel文件
#如有多个sheet可用“sheet”选项选择、“name”为输入名→"import"
#File/Url：显示Excel文件的存储路径
#Data Preview:用于预览Excel表，可以通过下拉菜单改变数据形式
#Import Options:Name：写R语言操作时，赋予Excel的新名字，如：example1
#Import Options:Sheet:如果Excel文件有多个Sheet表，在此选中需要操作的那个表，如：Sheet1
#Code Preview:这里自动生成操作相应的代码
#txt文件和csv文件的读取：
#方法1：代码
#txt文件的读取
example3<-read.table("文件名")
example3<-read.table("路径名")
example3
#csv文件的读取
example4<-read.csv("文件名")
example4<-read.csv("路径名")
example4
#方法二：鼠标操作
#最上方的"file"→"import dataset"→"from text(base)"→"browse"→选中操作文件→“Open”→“Import”
#“name”：写入输入名
#"heading"：指第一行是否作列名
#R语言文件的保存与读入
#保存代码1：saveRDS(R语言中文件名,"新文件名.保存文件格式")
#保存代码2：saveRDS(R语言中文件名,"~/文件夹1/.../文件夹n/新文件名.保存文件格式")
saveRDS(example1,"example1save1.rds")
#读取代码1：readRDS(R语言中文件名,"文件名.fileformat")
#读取代码2：readRDS(R语言中文件名,"~/文件夹1/.../文件夹n/文件名.fileformat")
newexample1<-readRDS("example1save1.rds")
#验证两个文件是否相同
identical(example1,newexample1)
##第二章-4.数据的概览和查看数据细节
#运行网络故障，可能是终端维护，可更换终端处理（视频的看不清楚，需要用再研究吧）
##查看数据的概览及细节
example1 <- read_excel("YKH.example1.xlsx")
View(example1)#查看表格的内容
View(example1[c(1:3),c(1:3)])#查看前3行、前3列的内容
class(example1)#查看数据是什么格式
class(example1[,1])#看第1列的数据是什么格式
class(example1[,3])#看第3列的数据是什么格式
dim(example1)#查看数据有几行、几列
nrow(example1)#查看数据有几行
ncol(example1)#查看数据有几列
head(example1)#查看数据的开头6行
tail(example1)#查看数据的最后6行
rownames(example1)[c(1:3)]#看第1-3行的行名
colnames(example1)[c(1:3)]#看第1-3列的列名
#第二章-5.查看缺失值、数据总结
#查看是否有缺失值(基础)：is.na() ; is.nan();
#对于矩阵或数据集查看某行/列是否有缺失值：an(is.na()) ;
#数据整体总结：summary() ; rowSums()colSums();
#进行额外计算：mutate()
#查看缺失值
#Na：表示缺失值，这里的值可以是数字，也可以是字符。“TRUE”表示有缺失值。
#NaN：表示缺失值，这里的值表示数字。TRUE表示有缺失值。
v1<-c(1:10)
is.na(v1)#查看v1向量中，有没有缺失数字或字符
is.nan(v1)#查看v1向量中，有没有缺失数字
any(is.na(example1[,1]))#查看第1列有没有缺失值
any(is.na(example1[,3]))#查看第1列有没有缺失值
any(is.na(example1[,4]))#查看第1列有没有缺失值
apply(example1, 1,function(x) any(is.na(x)))#查看整个文件中每一行是否有缺失值，1表示行
apply(example1, 2,function(x) any(is.na(x)))#查看整个文件中每一列是否有缺失值，2表示列
#去掉缺失值
example1.complete1<-example1[complete.cases(example1),]#去掉含有缺失值的行
#这个代码只能去除行，不能去除列，因为行是一个样本，一般都是样本有缺失值时去掉该样本
example1.complete1
dim(example1.complete1)#查看还剩几行几列
#查看数据的总结
#缺失值也会在总结中提示
summary(example1)
summary(example1.complete1)
#增加新的一列
#下面的“Tolal1”、“Tolal2”是赋予新的列的列名
install.packages("dplyr")
library("dplyr")
t1<-mutate(example1.complete1,Total1=Price*Kilogram)#定义这一列的数值为某两列相乘
t1
t2<-mutate(example1.complete1,Total2=Price^2)#定义这一列的数值为Price这一列的平方
t2
example1.sub<-example1.complete1[c(3:4)]#把文件中第2-4列的数据提取出来
View(example1.sub)
colSums(example1.sub)#把每一列相加(只能用于数字格式)
rowSums(example1.sub)#把每一行相加(只能用于数字格式)
##第二章-6.分类数据、连续数据
#分类数据
install.packages("readxl")
library("readxl")
example1 <- read_excel("YKH.example1.xlsx")
View(example1)
region1<-example1$Region#把文件“example1”中名称为“Region”这一列的内容提取出来
region1
region<-as.factor(region1)#把这一列数据由文本形式转化成分类变量形式
class(region)#可以看到，这一列的数据已经变成了分类变量形式：factor
levels(region)#查看这一列的唯一项（不能用于数字类型的数据）
summary(region)#分类变量的总结：查看这一列的各个唯一项分别有多少个
#连续数据（数字数据）
install.packages("readxl")
library("readxl")
example1 <- read_excel("YKH.example1.xlsx")
View(example1)
price<-example1$Price#把文件“example1”中名称为“Price”这一列的内容提取出来
class(price)#可以看到，这一列的数据为数字形式：numeric
sum(price)#计算这一列的总和
max(price)#找出这一列的最大值
min(price)#找出这一列的最小值
levels(price)#查看这一列的唯一项，发现失败了，计算结果显示Null，说明不能用于数字类型的数据
summary(price)#数字变量的总结：查看这一列的最大最小值、中位数、平均值、25%和75%的值
which.max(price)#寻找这一列数据中第几个是最大值
which.min(price)#寻找这一列数据中第几个是最小值
#如何把数字形式的数据转变为分类形式
kilo<-example1$Kilogram#把文件“example1”中名称为“Kilogram”这一列的内容提取出来
kilogram<-as.factor(kilo)#把数字形式数据转变为分类形式
class(kilogram)#可以看到，这一列的数据已经变成了分类变量形式：factor
levels(kilogram)#查看这一列的唯一项
summary(kilogram)#分类变量的总结：查看这一列的各个唯一项分别有多少个
##第二章-7.日期数据、时间数据
##定义
#① 日期数据(如："2018-09-19")
#② 时间数据(如："2018-09-19 20:00")
##系统日期数据符号
# %Y：4位年 eg：2007
# %y：2位年 eg：07
# %m：月 eg：01-12
# %B：月 eg：十月
# %b：月缩写 eg：10月
# %d：日 eg：01-31
# %A：星期几 eg：星期一
# %a：星期几缩写 eg：周一
##（1）日期数据的计算
install.packages("lubridate")#进行日期操作需要"lubridate"拓展包
library("lubridate")
Sys.Date()#显示系统日期（现在的日期）
format(Sys.Date(),format="%m %d %Y %a") #改变日期数据输出形式，"%m %d %Y %a"为指定的输出格式
#把文本数据转换成日期数据
date1<-as.Date("08/16/1975","%m/%d/%Y") #"08/16/1975"分别对应"%m/%d/%Y"
date1
#把数字数据转换成日期数据
as.Date(2053,origin="1970-01-01") #从1970-01-01开始2053天的日期，1970-01-01为sas日期数据的默认起
as.Date(-5252,origin="1990-01-01") #从1990-01-01之前2053天的日期，1990-01-01为excel日期数据的默认起点
##（2）时间数据的计算
#POSIXlt：把每一个时间数据保存成包括年份，月份，日期，时，分，秒等等元素的列表
#POSIXct：把时间数据保存成从起始点(1970-01-01 00:00.00)开始的秒数
#把时间字符串转换成时间数据
time1<-as.POSIXlt("2008-04-06 10:11:01",format="%Y-%m-%d %H:%M:%S")
time1#输出结果为"2008-04-06 10:11:01 HKT"，HKT表示时区（不做详述）
time2<-as.POSIXct(99999,origin="1970-01-01 00:00.00") #查看年月日时分秒等各项数据（月是0-11；年是从1900年开始数）
time2
#把年月日时分秒等各项数据作为列表展开并查看。
#unclass功能：月份显示3，是因为POSIXlt功能的月份是从0开始数的。
#unclass功能：年份显示108，是因为POSIXlt功能的年份是从1900开始数的。
unclass(time1)
time1$year#调取年分查看
time1$min#调取分钟查看
#更方便地从日期和时间数据当中提取所需要的成分的代码
install.packages("lubridate")
library("lubridate")
year(time1)#查看年
month(time1)#查看月
day(time1)#查看日
hour(time1)#查看时
minute(time1)#查看分
second(time1)#查看秒
wday(time1)#查看星期几，得出一个数字，1~7分别代表星期日~星期六
wday(time1,label=TRUE)#查看星期几
week(time1)#查看这是这一年的第几周
yday(time1)#查看是年中的第几天
mday(time1)#查看是月中的第几天
tz(time1)#查看时区
#更灵活地把文本转换成日期和时间数据的代码
#指定日期元素的顺序，可自动识别的分隔符包括:“-”、“/”、“.”以及无分隔符
#把各种形式的日期读取，并生成为如"2018-09-18"的格式
#ymd() 年，月，日
#ydm() 年，日，月
#mdy() 月，日，年
#dmy() 日，月，年
#hm() 时，分
#hms() 时，分，秒
#ymd_hms() 年，月，日，时，分，秒
ymd("2018-09-18")
dmy("18-09-2018")
mdy("09-18-2018")
ymd("2018/09/18")
ymd("2018\09\18")#“\”R语言识别不出，不能这样写，会报错
#生成一个包含姓名与日期的数据框
names<-c("Jack","Tom","Alice")#生成包含姓名的向量
DOB<-c("01/05/1990","08/16/2005","02/20/2000")#生成包含日期的向量
log_time<-c("2018-06-01 12:00.00","2017-01-09 08:20.30","2018-09-20 20:40.02")#生成包含时间的向量
dat<-data.frame(name=names,DOB=DOB,logtime=log_time)#把三个向量生成数据框
dat
str(dat)#查看结构，发现三列数据都是文本形式
dat$DOB<-as.Date(dat$DOB,"%m/%d/%Y")#所以，把名字为dat的数据框中的日期文本转变为日期数据
dat$logtime<-ymd_hms(dat$logtime)#所以，把名字为dat的数据框中的时间文本转变为时间数据
dat
str(dat)#再查看结构，发现相应数据的格式已经转变为日期数据格式和时间数据格式
dat$logy-year(dat$logtime)#提取年份
dat$logm <- month(dat$logtime)#提取月份
dat$logwdayy<- wday(dat$logtime, label=TRUE)#提取星期几
dat#查看，与上面3个代码一起run
##第二章-8.如何处理数据中的文本
#创建一个文本
as.character("Hello,world")#文本内容用双引号或者单引号括起来都可以
as.character('Hello,world')
#判断是否为文本
is.character("Hello,world")
is.character(1)
#判断文本数据的长度
nchar("Hello")
nchar("world")
nchar("Hello,world")#逗号也算一个文本
#大写变小写
tolower("HELLO,WORLD")
#小写变大写
toupper("hello,world")
#文本的拼接（不写分隔符，默认为空格）
paste("Hello","world",sep=" ")#sep后面引号内的内容，为连接两个文本之间的分隔符
paste("Hello","world",sep=",")
paste("Hello","world",sep="傻逼")
paste("Hello","world")#没设定分隔符，两个文本之间的分隔符默认为空格
#另一种拼接方法 paste0 （无缝拼接）
paste0("Hello","world")
paste("Hello","world",sep="")#上面的方法也能做到无缝拼接
#文本的剪切
substr("Hello,world",1,5)#“1,5”为第一个剪切的文本的位置、最后一个剪切文本的位置
substr("Hello,world",7,12)
#另一种剪切方式（分割）
strsplit("Hello,world",split=",")#把某个部分","去掉
strsplit("Hello,world,",split=",")#可以看到去掉了所有的","
#文本的替换(只能替换相同长度的文本)
x<-c("Hello,,world")
substr(x,6,7)<-":)" #把,,替换成了:)
x
x1<-c("Hello,world")
substr(x1,7,11)<-"Beijing"
x1 #结果可以看到，只把world替换成了Beijing的前5个字符
#想要把不同长度的文本相互替换，只能先把文本剪切成几个不同的部分，再粘合到一起
x2<-c("Hello,world")
x3<-substr("Hello,world",1,6)
x4<-c("Beijing")
paste0(x3,x4)
#规律替换(只替换第1个)
sub("A","a","AAaaBbBb")
sub("b","B","AAaaBbBb")
#规律替换(全部替换)
gsub("A","a","AAaaBbBb")
gsub("b","B","AAaaBbBb")
##第二章-9.宽数据与长数据的转换
#①数据结构：包括宽数据（wide fotmat）、长数据(long format)
#②数据结构的重组 (数据重组/reshape)：
#安装加载reshape2、data.table两个R包
install.packages("reshape2")
install.packages("data.table")
library("reshape2")
library("data.table")
load("diagnosis.RData")#打开操作文件
str(data.wide)#查看宽数据有什么变量（可在右上角操作）
View(data.wide)#查看宽数据表格
##把宽数据转换成长数据
#data：要转换的宽数据
#id.var：保持不变（保留）的数据
#measure.var：要转换的数据
#variable.name：转换后的分类变量名
#value.name：转换后的数值变量名
data.long <- reshape2::melt(data = data.wide,
id.var = c('ID','Age','Sex','Status'),
measure.var=c('TestA', 'TestB', 'TestC'),
variable.name='Test',
value.name='Test.Result')
View(data.long)
data.wide[data.wide$ID==1,]#查看转换前的宽数据的第1个人的数据
data.long[data.long$ID==1,]#查看转换后的长数据的第1个人的数据
##把长数据转换成宽数据
#data：要转换的长数据
#formula：~前为保持不变的变量；~后为分类变量的转换依据
#value.var：转换后分类变量的值的来源
data.wide1<-reshape2::dcast(data=data.long,
formula=ID+Age+Sex+Status~Test,
value.var = 'Test.Result')
View(data.wide1)
data.wide1[data.wide$ID==1,]#查看转换后的第1个人的数据
##重组多个数据：
#因为分类变量相关的变量不仅有数值，还有日期，如果多个分量变量相关数据都保留，该如何操作
#(下面三段方法不想看，要用的时候再说吧)
#数据重组方法1：根据分类数据，选择相应的列
data.long.multi <- reshape2::melt(data = data.wide,
id.vars=c('ID','Age','Sex','Status','TestA.date','TestB.date',
'TestC.date'),
measure.vars = c('TestA', 'TestB','TestC'),
variable.name='Test',
value.name='Test.Result')
data.long.multi$Test.date<-ifelse(data.long.multi$Test=='TestA',data.long.multi$TestA.date,
ifelse(data.long.multi$Test=='TestB',data.long.multi$TestB.date,data.long.multi$TestC.date))
# tell R Test.date is a date variable
data.long.multi$Test.date<-as.Date(data.long.multi$Test.date,origin='1970-01-01')
# drop TestA.date, TestB.date, TestC.date
data.long.multi <-
subset(data.long.multi,select=c(TestA.date,TestB.date,TestC.date))
head(data.long.multi)
#数据重组方法2：自己动手，先剪切，再拼接
# A
data.wide.A<- subset(data.wide,select=c('ID','Age','Sex','Status','TestA','TestA.date'))
names(data.wide.A)[names(data.wide.A) == 'TestA'] <- 'Test.Result'
names(data.wide.A)[names(data.wide.A) == 'TestA.date'] <-'Test.date'
data.wide.A$Test <- 'TestA'
# B
data.wide.B<- subset(data.wide,select=c('ID','Age','Sex','Status','TestB','TestB.date'))
names(data.wide.B)[names(data.wide.B) == 'TestB'] <- 'Test.Result'
names(data.wide.B)[names(data.wide.B) == 'TestB.date'] <- 'Test.date'
data.wide.B$Test <- 'TestB'
# C
data.wide.C<- subset(data.wide,select=c('ID','Age','Sex','Status','TestC','TestC.date'))
names(data.wide.C)[names(data.wide.C) == 'TestC'] <- 'Test.Result'
names(data.wide.C)[names(data.wide.C) == 'TestC.date'] <-'Test.date'
data.wide.C$Test <- 'TestC'
# combine A,B,C
data.long.multi2 <- rbind(data.wide.A,data.wide.B,data.wide.C)
head(data.long.multi2)
#数据重组方法3：用data.table包(简单点)
# data.wide to data.table format
data.wide.table <- as.data.table(data.wide)
str(data.wide.table)
View(data.wide.table)
data.long.multi3 <- data.table::melt(data = data.wide.table,
id.vars = c('ID','Age','Sex','Status'),
measure.vars = list(c('TestA', 'TestB','TestC'), c('TestA.date','TestB.date','TestC.date')),
variable.name = 'Test',
value.name = c('Test.Result', 'Test.date'))
View(data.long.multi3)
head(data.long.multi3)
###第三章：数据分析和绘图
##第三章-1.常用统计方法和描述统计
#（1）分析一个变量的分布：①分类变量；②连续变量
#（2）分析两个变量的关系：①分类变量&分类变量；②连续变量&分类变量；③连续变量&连续变量
#（3）描述统计和基本特征表：
#安装并加载需要的安装包
install.packages('table1')
install.packages('lubridate')
library(table1)
library(lubridate)
#加载需要使用的文件
load('diagnosis.RData')
View(data.wide)
#给每一个变量加上label(说清楚每一个数值对应的含义)
data.wide$Status <- factor(data.wide$Status,
levels=c(0,1),
labels=c("No Disease", "Disease"))
data.wide$Sex <- factor(data.wide$Sex,
levels=c(0,1),
labels=c("Male", "Female"))
#attach可以使接下来的操作都是用data.wide数据，不用每次都重复输入了
attach(data.wide)
### 描述变量的分布
#描述分类变量的分布（例如：性别(男、女)、健康状况(患病、健康)……）
#描述分类变量分布的常用的统计量：频数、百分比
##操作如下：
#描述频数
mytable_Sex <- table(Sex)
mytable_Sex
#描述百分比
prop.table(mytable_Sex)
#百分比假设检验
prop.test(mytable_Sex,p=0.5)#假设某分类变量百分比为50%，检验假设是否成立
#描述连续变量的分布（例如：年龄、检查指标值……）
#描述连续变量分布的常用的统计量：均值、标准差，中位数四分位距、极差
#当变量服从正态分布时使用：t检验
#当变量不服从正态分布时使用：Wilcoxon符号秩检验
##操作如下：
#看连续变量的概括
summary(Age)#包含：最小值、25%分位数、中位数、75%分位数、最大值
#检查某个均值是否为假设的常数（t检验：用于变量服从正态分布）
t.test(Age,mu=35)
#检查某个均值是否为假设的常数（Wilcoxon符号秩检验：用于变量不服从正态分布）
wilcox.test(Age,mu=35)
##组间比较
#（1）分类变量&分类变量（例如:患病和性别的关系）
#分析方法:列联表和卡方检验(x’test )
#列联表
mytable <- table(Sex,Status)
mytable
#卡方检验（eg：检验性别对患病情况是否有显著影响；看p值）
chisq.test(mytable)
#统计每一行变量的情况
margin.table(mytable, 1)
#统计每一列变量的情况
margin.table(mytable, 2)
#(2)连续变量（正态分布）&分类变量（二分类）(例如:年龄和性别的关系)
#分析方法:t检验
#假设方差相等时，查看t检验统计量（t）、p值、两组分别的均值
t.test(Age ~ Sex,var.equal = TRUE)
#假设方差不相等时，查看t检验统计量（t）、p值、两组分别的均值
t.test(Age ~ Sex,var.equal = FALSE)
#(3)连续变量（非正态分布）&分类变量（二分类）(例如:检验的指标值和性别的关系)
#分析方法:Mann-Whitney U检验(Wilcoxon秩和检验)
#用Mann-Whitney U检验（Wilcoxon秩和检验）（比较的是两组的分布是否一致，而不是均值；结果看p值）
wilcox.test(TestA ~ Sex)
#(4)连续变量（正态分布）&分类变量（多分类）(例如:检查指标值(log转换后)和检查年度的关系)
#分析方法:单变量方差分析(ANOVA)
#把检查年度生成多分类变量
data.wide$TestA.year<-as.factor(ifelse(year(data.wide$TestA.date)<2009,"<2009",
ifelse(year(data.wide$TestA.date)==2009,"2009",">2009")))
#用单变量方差分析（ANOVA）进行分析，查看样本均值、每组均值及频数、F检验统计量和p值
#log(TestA) ：可以使连续变量更符合正态分布
fit <- aov(log(TestA) ~ TestA.year)
model.tables(fit,"means")#查看样本均值、每组均值及频数
summary(fit)#查看F检验统计量、p值
#(5)连续变量（非正态分布）&分类变量（多分类）(例如:检查指标值和检查年度的关系)
#分析方法:Kruskal-Wallis检验（结果看p值）
kruskal.test(TestA ~ TestA.year)
#(6)连续变量&连续变量(例如:年龄，检查指标值)
#分析方法:①均为正态分布：Pearson相关系数；②不均为正态分布：Spearman相关系数
#均符合正态分布，用Pearson相关系数分析
#查看Pearson相关系数
cor(TestA,TestB, method="pearson")
#若还需看p值（一般不需要）
cor.test(TestA,TestB, method="pearson")
#不均符合正态分布，用Spearman相关系数
#查看Spearman相关系数
cor(TestA,TestB, method="spearman")
##描述统计和基本特征图
#需要的扩展包：table1
library(table1)
#用table1包，自动生成基本特征表
table1(~ Status + Sex + Age + TestA + TestB + TestC + TestA.year, data=data.wide)
#把患病情况分组展示数据（把分类变量放竖线后面）
table1(~ Sex + Age + TestA + TestB + TestC + TestA.year| Status, data=data.wide)
#在table 1增加更多的数据
table1(~ Sex + Age + TestA + TestB + TestC + TestA.year| Status, data=data.wide,
render.continuous=c(.="Mean (SD)", .="Median [Min, Max]",.="Median [Q1, Q3]",.="Median [IQR]"))
#在table 1增加p值方法如下：
load('diagnosis.RData')#加载表格
#把状态、性别定义为分类变量
data.wide$Status <- factor(data.wide$Status,
levels=c(0, 1, 2),
labels=c("No Disease", "Disease", "P-value"))
data.wide$Sex <- factor(data.wide$Sex,
levels=c(0,1),
labels=c("Male", "Female"))
#显示组间比较的p值
data <- data.wide
outcome <- data.wide$Status
rndr <- function(x, name, ...) {
if (length(x) == 0) {
y <- data[[name]]
s <- rep("", length(render.default(x=y, name=name, ...)))
if (is.numeric(y)) {
p <- t.test(y ~ outcome)$p.value
} else {
p <- chisq.test(table(y, droplevels(outcome)))$p.value
}
s[2] <- sub("<", "<", format.pval(p, digits=3, eps=0.001))
s
} else {
render.default(x=x, name=name, ...)
}
}
rndr.strat <- function(label, n, ...) {
ifelse(n==0, label, render.strat.default(label, n, ...))
}
#生成含有p值的表格
table1(~ Sex + Age + TestA + TestB + TestC | Status, data=data.wide,
render.continuous=c(.="Mean (SD)"),droplevels=F, render=rndr, render.strat=rndr.strat, overall=F)
#这个意思是清除前面所有操作
rm(list = ls())
##第三章-2.基本统计图
#频数:条形图
#百分比:饼图
#连续变量的分布:直方图、箱线图
#分组比较:多组箱线图
#相关系数:散点图
#这个意思是清除前面所有操作
rm(list = ls())
#条形图（分类变量的频数）
load('diagnosis.RData')
gender<-table(data.wide$Sex)
glables<-c('男人','女人')#定义分类变量名称
barplot(gender)#画出条形图
barplot(gender,names.arg= glables)#在条形图加上分类变量名称
barplot(gender,names.arg= glables,main='性别人数对比图')#在条形图加上标题
barplot(gender,names.arg= glables,main='性别人数对比图',xlab='sex',ylab='count')#在条形图加上x、y坐标轴名称
#扇形图
gender<-table(data.wide$Sex)
glables<-c('男人','女人')
pie(gender)#画出扇形图
pie(gender,glables, main='性别人数对比图')
#直方图（连续变量）
hist(data.wide$Age)
hist(data.wide$Age,main='性别人数对比图',xlab='年龄')
#“breaks-5”把显示的坐标轴（20-70）上的直方图分成5条柱子
hist(data.wide$Age,main='性别人数对比图',xlab='年龄',breaks=5)
#规定直方图在“15,25,35,45,55,65,75”这几个点分栏
hist(data.wide$Age,main='性别人数对比图',xlab='年龄',breaks=c(15,25,35,45,55,65,75))
#箱线图（连续变量）
boxplot(data.wide$Age)
boxplot(data.wide$Age,main='性别人数对比图',ylab='年龄')
#画水平参考线（高度60，颜色红色）
abline(h=60,col='red')
#把箱线图转成打横（连续变量）
boxplot(data.wide$Age,horizontal=TRUE,main='性别人数对比图',ylab='年龄')
#画垂直参考线（高度60，颜色红色）
abline(v=60,col='blue')
#分组箱线图（连续变量）
boxplot(Age~Sex,data=data.wide)#~的左边是用来画箱线图的变量，右边是分组的变量
boxplot(Age~Sex,data=data.wide,names=glables,main='性别人数对比图',xlab = '性别',ylab = '年龄')
#散点图
plot(data.wide$Age,data.wide$TestC)
plot(data.wide$Age,data.wide$TestC,main='各年龄段TestC评分',xlab= '年龄',ylab = 'TestC评分')
#按照y=a+bx函数给a、b赋值画出参考线
abline(a=-20,b=1,col='green')
#在特定位置加上文本
text(x=40,y=25,'y=-20+x',col="green")
#在1个图片中同时画多个图
par(mfrow=c(2,2))#把图片分成2行、2列的4个区域
#这个时候再画图，图片就会分别放在4个区域内
barplot(gender,names.arg= glables,main='性别人数对比图',xlab='性别',ylab='人数')
pie(gender,glables, main='性别人数对比图')
boxplot(data.wide$Age,main='性别人数对比图',ylab='年龄')
abline(h=60,col='red')
plot(data.wide$Age,data.wide$TestC,main='各年龄段TestC评分',xlab= '年龄',ylab = 'TestC评分')
abline(a=-20,b=1,col='green')
text(x=40,y=25,'y=-20+x',col="green")
#保存图片
#方法1（鼠标操作）：右下角框内→“Export”→“Save as Image”→调整参数→“Save”
#调整参数：lmage format：更改图片格式
#调整参数：Directory：更改保存位置
#调整参数：File name：设置图片保存名称
#调整参数：Width、Height：调整长宽（长宽不够，图片会显示不全）
#调整参数：Update Preview：更新图形并预览
#方法2（代码）
#保存图片常用代码：png(),jpeg(),bmp(),tiff()
#以下以保存为png格式图片演示
#图片名必须是“名字.格式”这种形式
png(filename = '图形.png',width =3000,height = 1500,res=300)#res是像素
par(mfrow=c(2,2))
barplot(gender,names.arg= glables,main='性别人数对比图',xlab='性别',ylab='人数')
pie(gender, glables,main='性别人数对比图')
boxplot(data.wide$Age,main='性别人数对比图',ylab='年龄')
abline(h=60,col='red')
plot(data.wide$Age,data.wide$TestC,main='各年龄段TestC评分',xlab= '年龄',ylab = 'TestC评分')
abline(a=-20,b=1,col='green')
text(x=40,y=25,'y=-20+x',col="green")
dev.off()#保存：运行结果显示“RStudioGD=2”表示已正常保存
##第三章-3.利用ggplot2绘制漂亮的图
#清除之前的痕迹
rm(list = ls())
#安装加载相应的R包
install.packages('colorspace')
library(colorspace)
install.packages('ggplot2')
library(ggplot2)
#加载数据
load('diagnosis.Rdata')
View(data.wide)
#把二分类变量的1、0转换成文字
data.wide$Status<- factor(data.wide$Status,
levels=c(0,1),
labels=c("没病", "有病"))
data.wide$Sex <-factor(data.wide$Sex,
levels=c(0,1),
labels=c("Male", "Female"))
#①条形图（分类变量的频数）
theme_update(plot.title=element_text(hjust= 0.5))#设定题目放在中间
qplot(Sex,data=data.wide,geom='bar' )#data=data.wide表示用哪个文件；geom='bar'表示制作的格式是条形图；
qplot(Sex,data=data.wide,geom='bar',main='性别条形图' )
qplot(Sex,data=data.wide,geom='bar',main='性别条形图' ,ylab='频数')
#②直方图
qplot(Age,data =data.wide,geom='histogram' )#第一个数据“Age”为横坐标，轴坐标为频数；
#调整每一格矩形的宽度
qplot(Age,data =data.wide,geom='histogram',binwidth=10 )
qplot(Age,data =data.wide,geom='histogram',binwidth=5 )
qplot(Age,data =data.wide,geom='histogram',binwidth=1 )
#写出图片名、x轴、y轴的名称
qplot(Age,data =data.wide,geom='histogram',binwidth=5,
main='年龄的直方图',xlab = '年龄（年）',ylab = '频数' )
#设置x、y轴的长度
qplot(Age,data =data.wide,geom='histogram',binwidth=5,
main='年龄的直方图',xlab = '年龄（年）',
ylab = '频数',xlim = c(10,80) ,ylim = c(0,50))
#改变填充颜色
qplot(Age,data =data.wide,geom='histogram',binwidth=5,
main='年龄的直方图',xlab = '年龄（年）',
ylab = '频数',xlim = c(10,80) ,ylim = c(0,50),fill=I('blue'))#必须要加“I”，软件才能识别出是颜色指令
#改变线条颜色
qplot(Age,data =data.wide,geom='histogram',binwidth=5,
main='年龄的直方图',xlab = '年龄（年）',
ylab = '频数',xlim = c(10,80) ,ylim = c(0,50),fill=I('blue'),col=I('red'))
#改变颜色透明度
qplot(Age,data =data.wide,geom='histogram',binwidth=5,
main='年龄的直方图',xlab = '年龄（年）',
ylab = '频数',xlim = c(10,80) ,ylim = c(0,50),
fill=I('blue'),col=I('red'),alpha=I(0.5))
#③箱线图（连续变量）
#画出一个不分组的箱线图
qplot(x='1',y=Age,data=data.wide,geom='boxplot' )
qplot(x='1',y=Age,data=data.wide,geom='boxplot',xlab = 'Age' )#不定义x轴的话，x轴下面写1，不好看；
#画分组的箱线图
qplot(x=Sex,y=Age,data=data.wide,geom='boxplot' )
#按分组填充不同颜色
qplot(x=Sex,y=Age,data=data.wide,geom='boxplot',fill=Sex)
#去掉上下的点（异常值？）
qplot(x=Sex,y=Age,data=data.wide,geom='boxplot',fill=Sex,outlier.shape=NA)
#加上散点图
qplot(x=Sex,y=Age,data=data.wide,geom=c('boxplot','jitter'),fill=Sex,outlier.shape=NA)
#④做分组箱线图：
#加载扩展包
install.packages(reshape2)
library(reshape2)
#把宽数据转成长数据
data.long <-reshape2::melt(data = data.wide,
id.vars =c('ID','Age','Sex','Status'),
measure.vars =c('TestA', 'TestB','TestC'),
variable.name='Test',
value.name='Test.Result')
#作图
qplot(x=Test,y=Test.Result, data=data.long, geom="boxplot", fill=Test)
#⑤散点图：
qplot(x=Age,y=TestC,data=data.wide,geom='point' )
#标出x、y轴的名称
qplot(x=Age,y=TestC,data=data.wide,geom='point',ylab = 'TestC的分数',xlab = '年龄' )
#调节点的大小、形状
qplot(x=Age,y=TestC,data=data.wide,geom='point',ylab = 'TestC的分数',xlab = '年龄',
size=I(3),shape=I('triangle') )
#给不同组的点设置不同的颜色和形状
qplot(x=Age,y=TestC,data=data.wide,geom='point',ylab = 'TestC的分数',xlab = '年龄',
size=I(3) ,col=Sex,shape=Sex)
#⑥把一个图分开几行几列表示
#把散点图按性别分行，按患病情况分列表示
qplot(x=Age,y=TestC,data=data.wide,geom='point',ylab = 'TestC的分数',xlab = '年龄',
facets =Sex~Status)
#调整点的颜色、形状
qplot(x=Age,y=TestC,data=data.wide,geom='point',ylab = 'TestC的分数',xlab = '年龄',
facets =Sex~Status,col=Sex,shape=Sex)
#只按性别分行
qplot(x=Age,y=TestC,data=data.wide,geom='point',ylab = 'TestC的分数',xlab = '年龄',
facets =Sex~.)
#只按患病情况分列
qplot(x=Age,y=TestC,data=data.wide,geom='point',ylab = 'TestC的分数',xlab = '年龄',
facets =.~Status)
#其他图如箱线图分组展示也是同样的操作
qplot(x=Sex,y=Age,data =data.wide,geom='boxplot',
facets =Sex~Status )
#⑦ ggplot源画图比上面的扩展包更具优势：①删减内容方便；②可以进行精细图片的调整；
#作图
ggplot(aes(x=Sex,y=Age),data=data.wide)+
geom_boxplot(aes(fill=Sex))+
facet_grid(.~Status)+
ggtitle('根据性别和患病情况的年龄分组')
#加上散点图：geom_jitter(position =position_jitter(0.4),size=1；
#加上y轴命名及范围限制：scale_y_continuous('年龄（年）',limits = c(10,80)；
ggplot(aes(x=Sex,y=Age),data=data.wide)+
geom_boxplot(aes(fill=Sex),outlier.shape =NA)+
geom_jitter(position =position_jitter(0.4),size=1)+
facet_grid(.~Status)+
ggtitle('根据性别和患病情况的年龄分组')+
scale_y_continuous('年龄（年）',limits = c(10,80))
##第三章-4.线性回归模型
##模型形式:因变量Y=(多个)自变量X的线性组合+残差项Ei
##模型形式:Y=β+β1*Xi1+β₂*Xi2+...+ βp,Xip+Ei
##适用范围:
#①因变量Y是连续变量，if not →第19课，广义线性模型
#②每行数据之间相互独立，如：时间序列，重复测量，纵向数据等适不适合使用线性回归模型的
#③残差分布独立于自变量X →第16课，模型诊断和改进
#④因变量Y和自变量X之间是线性关系 →第17课，样条函数
#准备工作
rm(list = ls())
install.packages('ggplot2')
library(ggplot2)
load('diagnosis.RData')
View(data.wide)
data.wide$Status <- factor(data.wide$Status,
levels=c(0,1),
labels=c("没病", "有病"))
data.wide$Sex <- factor(data.wide$Sex,
levels=c(0,1),
labels=c("男的", "女的"))
attach(data.wide)#以下操作都默认用data.wide数据
# 如果没有进行“attach(data.wide)”这一步，编写的代码应该如下，
# model1 <- lm(TestC ~ Age, data=data.wide)
##线性回归
# 1.1 一个自变量
model1 <- lm(TestC ~ Age)#~之前为因变量，之后为自变量
model1
# 1.2 多个自变量
model2 <- lm(TestC ~ Age + TestA + TestB)
model2
# 1.3 分类自变量 (需要引入哑变量)
model3 <- lm(TestC ~ Age + relevel(Sex,"Male")) #以男性作为参考组
model3
model3.1 <- lm(TestC ~ Age + relevel(Sex,"Female")) #以女性作为参考组
model3.1
# 1.4 交互项
#交互项：A*B = A + B + A:B
# "*" ：会输出两个自变量及两个自变量的交互项的回归系数
model4 <- lm(TestC ~ Age*relevel(Sex,"Male"))
model4
# ":" ：只输出两个自变量的交互项的回归系数
model4.1 <- lm(TestC ~ Age + Age:relevel(Sex,"Male"))
model4.1
##查看模型表现需要看的值：F检验、模型整体p值、R平方、调整后R平方（运行结果的最后两行）
summary(model4) # F test p value, R^2, adj R^2
##如果有需要，还可以通过AIC和BIC功能查看AIC值和BIC值
AIC(model4)
BIC(model4)
#回归系数
#查看模型的：回归系数、标准误、p值
summary(model4)$coefficients
#查看模型的回归系数的95%置信区间
confint(model4, level=0.95)
#把模型的回归系数提取出来（包括截距、各个自变量的回归系数 ）
coef(model4)
#查看4个回归系数的第2个
coef(model4)[2]
# 3. 用回归模型做预测
#将得出的截距和回归系数放进方程可得到公式
#将新的自变量X的数据带入回归模型的公式中，得到因变量的预测值
# 用得出回归系数的患者数据放进得出的公式中，求得因变量的预测值
predict(model4)
#用新的（不是得出回归系数的）患者数据放进公式中，求得因变量的预测值
#创建新的患者的数据集
new.patient <- data.frame(Age=c(30, 35, 40),Sex=c(0,1,1))
new.patient$Sex <- factor(new.patient$Sex,
levels=c(0,1),
labels=c("Male", "Female"))
new.patient
#用新的（不是得出回归系数的）患者数据放进公式中，求得因变量的预测值
new.patient$TestC_hat <- predict(model4,newdata = new.patient)#$TestC_hat：这是new.patient数据框中的一个新列，用来存储预测值
new.patient
#用回归模型的预测结果的图，展示因变量与自变量间的关系
#画出散点图
#coef(model1)[1]：model1的截距
#coef(model1)[2]：model1的回归系数（斜率）
plot_model1 <- qplot(x=Age, y=TestC, data=data.wide, geom = "point")+
geom_abline(intercept = coef(model1)[1], slope = coef(model1)[2])
plot_model1
#把男女的散点分开，并做出各自的线
#coef(model3)[3]：Sex（非"Male"组）的系数，表示与"Male"组相比，其他组在TestC上的平均差异，所以在截距部分要加上coef(model3)[3]
model3 <- lm(TestC ~ Age + relevel(Sex,"Male"))
model3
plot_model3 <- qplot(x=Age, y=TestC, data=data.wide, geom = "point", col=Sex)+
geom_abline(intercept = coef(model3)[1], slope = coef(model3)[2], col="red") +
geom_abline(intercept = coef(model3)[1]+coef(model3)[3], slope = coef(model3)[2], col="blue")
plot_model3
#把男女的散点分开，并做出各自的线
#coef(model4)[4]：Age 和 Sex 之间的交互项系数，表示在非 Male 组中，Age 每增加一个单位，TestC 相比 Male 组增加的额外数量
model4 <- lm(TestC ~ Age*relevel(Sex,"Male"))
model4
plot_model4 <- qplot(x=Age, y=TestC, data=data.wide, geom = "point", col=Sex)+
geom_abline(intercept = coef(model4)[1], slope = coef(model4)[2], col="red") +
geom_abline(intercept = coef(model4)[1]+coef(model4)[3], slope = coef(model4)[2]+coef(model4)[4], col="blue")
plot_model4
###第三章-5.回归模型的诊断和改进
#操作前准备
rm(list = ls())
install.packages('car')
install.packages('ggplot2')
library(car)
library(ggplot2)
load('diagnosis.RData')
View（data.wide）
data.wide$Status <- factor(data.wide$Status,
levels=c(0,1),
labels=c("No Disease", "Disease"))
data.wide$Sex <- factor(data.wide$Sex,
levels=c(0,1),
labels=c("Male", "Female"))
attach(data.wide)
## 1. 检验数据的正态性
# 1.1 统计检验（Shapiro-Wilk test）
#原假设：H0: 数据来自于正态分布总体。
#但是正态性受样本量影响比较大，样本量小的时候，p值偏大，不容易拒绝原假设。
#而样本量大的时候，又太敏感，太容易拒绝原假设。
#所以，不推荐直接根据检验统计量及p值判断数据是否服从正态分布。
#可以推荐图形检验（直方图，Q-Q图）判断数据是否服从正态分布。
shapiro.test(TestC)
shapiro.test(log(TestC))
# 1.2 带有密度曲线的直方图
#直方图如果对称，考虑服从正态分布。
hist(TestC,probability=T)#画直方图
lines(density(TestC),col=2)#加上密度曲线
#有些不符合正态分布的数据，转换成对数函数后，变得符合正态分布了
hist(log(TestC),probability=T)
lines(density(log(TestC)),col=2)
#1.3 Q-Q图
#Q-Q图的点越接近正态分布参考线，越考虑服从正态分布。
qqnorm(TestC)
qqline(TestC)
#有些不符合正态分布的数据，转换成对数函数后，变得符合正态分布了
qqnorm(log(TestC))
qqline(log(TestC))
## 2. 回归模型的判断
# 2.1 评价模型表现指标：R平方（回归模型的决定系数）、调整后R平方
#R平方：因变量的方差中，可由自变量解释的部分所占的比例。
#调整后的R平方：因为放进回归模型的因变量越多，R平方的数值越大。所以不能仅仅通过R平方判断，纳入因变量的数量进行分析后，得出调整后的R平方。
#不同学科、不同的研究，对R平方的要求是不一样的。
#我们只能说，R平方越小，模型的改进空间越大，但在实际应用中，不需要追求R平方达到特定的数值。
#以下为求R平方、调整后R平方的代码
model2 <- lm(TestC ~ Age + TestA + TestB)
summary(model2)
# 2.2 线性关系
#红色的线越接近水平线，说明残差越不容易随预测值的变化而变化，表明具有较好的线性关系
plot(model2,1)
# 2.3 残差正态性
#观察残差与正态分布参考线的吻合程度，越吻合说明越符合状态分布
plot(model2,2)
# 2.4 方差齐性
#红色线不应有明显的趋势。红色的线越接近水平线，说明具有较好的线性关系
plot(model2,3)
# 2.5 强影响点
#指标：Cook距离
#Cook距离最大的前三个点会被标注
#Cook距离>1的点，可以认为是强影响点
plot(model2,4)#查看有影响的观测值图
plot(model2,5)#查看Cook's距离图，用于识别具有较大影响的观测值
# 2.6 多重共线性
#指标：用方差膨胀因子（VIF）评估，VIF越大，多重共线性越严重。
#一般认为：方差膨胀因子（VIF）小于10，不存在多重共线性。
vif(model2)
##模型的改进
#R平方值偏低 → 寻找被遗漏的自变量
#残差正态性 → 对不服从正态分布的变量进行转换（如前面提过的取对数（log））
#方差齐性 → 使用加权最小二乘法解决
#强影响点 → 剔除强影响点，并比较模型结果
#多重共线性 → 去除与其他自变量强相关的自变量
# 3. 请尝试其他方法来改进模型
#反正就是用下面两个回归模型分别去做，不具体演示了
model2_1 <- lm(TestC ~ Age + TestA + TestB + relevel(Sex,"Male"))
model2_2 <- lm(log(TestC) ~ Age + log(TestA) + log(TestB))
###第三章-6.回归模型的比较和变量筛选
#操作前准备
#去除之前的操作痕迹
rm(list = ls())
## install packages
install.packages('ggplot2')
install.packages('leaps')
install.packages('car')
install.packages('glmnet')
install.packages('QuantPsyc')
install.packages('relaimpo')
library(ggplot2)
library(leaps)
library(car)
library(glmnet)
library(QuantPsyc)
library(relaimpo)
## load data
load('diagnosis.RData')
## add labels
data.wide$Status <- factor(data.wide$Status,
levels=c(0,1),
labels=c("No Disease", "Disease"))
data.wide$Sex <- factor(data.wide$Sex,
levels=c(0,1),
labels=c("Male", "Female"))
View(data.wide)
## generate some new variabels
data.wide$TimeInterval_AC <- as.numeric(data.wide$TestC.date-data.wide$TestA.date)
data.wide$TimeInterval_BC <- as.numeric(data.wide$TestC.date-data.wide$TestB.date)
# attach data.wide
attach(data.wide)
## 1. 比较两个回归模型
# 1.1 模型是嵌套关系（一个模型完全在另一个模型中：Model2 = Model1 + X）
# anova命令：是前面介绍过的单变量方差分析同样的命令
model1 <- lm(TestC ~ Age)
model2 <- lm(TestC ~ Age + TestA + TestB)
anova(model1,model2)#主要查看F值；其次查看p值（若p值＜0.05,则说明加入TestA、TestB对模型有改善）
# anova(model1,model2, test = "Chisq") # can also use chi-squred test
# 1.2 模型不是嵌套关系（Model2 = Model1 + X1；Model3=Model1 + X2）
model2 <- lm(TestC ~ Age + TestA + TestB)
model3 <- lm(TestC ~ Age + relevel(Sex,"Male"))
AIC(model2,model3)#AIC越小，模型越好
## 2. 模型变量的选择
# 2.1 逐步回归（stepwise selection）
#基本思想：每次加入一个新的变量，或者去掉一个模型中现有的变量，较之前的模型是否改善，再决定是否加入或者剔除变量
#逐步回归简单易用，但是有很多缺陷，现在已经基本不推荐逐步回归法了
#需要使用的指令：stepAIC()
#建立一个完整的模型，这个模型需要包含所有的候选变量
model.full <- lm(TestC ~ Age + relevel(Sex,"Male") + relevel(Status,"No Disease") +
TestA + TestB + TimeInterval_AC + TimeInterval_BC)
#向前逐步回归法
model.forward <- stepAIC(model.full, direction="forward")#运行结果可以看到筛选的过程
summary(model.forward)#查看最后有几个变量剩下
#向后逐步回归法
model.backward <- stepAIC(model.full, direction="backward")
summary(model.backward)
#双向逐步回归法
model.both <- stepAIC(model.full, direction="both")
summary(model.both)
# 2.2 全部子集回归（all subsets）
# p是原始完整模型中变量个数。每个变量只有纳入或排除两种可能。
# 所以p个变量，共有（𝒑的二次方）个回归模型，去除所有变量都不纳入的回归模型后，还剩下（𝒑的二次方 − 𝟏）个模型。
# 使用的命令：regsubsets(y~x，nbest=k )，k是指：对于每种变量个数的模型，展示最好的k个模型
# 全部子集回归法可以根据预先设定的最终模型中变量数选择模型，或者选择拟合表现最佳的模型作为最终模型
model.all <- regsubsets(TestC ~ Age + relevel(Sex,"Male") + relevel(Status,"No Disease") +
TestA + TestB + TimeInterval_AC + TimeInterval_BC, nbest=3,data=data.wide)
# 得出的图是把最后得出的所有模型，按照模型的调整后的R平方值排序展示
# scale="adjr2"：根据调整后的R平方排序
# 按右下角框内的zoom可以放大图片查看
plot(model.all, scale="adjr2")
# 2.3 Lasso回归法（这个方法高级）
# Lasso: Least absolute shrinkage and selection operator（通过缩小回归系数，达到变量选择的目的）
# 中心思想：回归系数β越小越好
# 目标：使残差平方的均值加λ（Lambda）倍回归系数绝对值的和最小化的β
# 优势：可以避免过度拟合的问题，同时通过将一些β缩减为0，达到变量筛选的目的
# 把分类变量转换为数值变量（Lasso回归法认不出分类变量）
data.wide$Sex_num <- as.numeric(data.wide$Sex)
data.wide$Status_num <- as.numeric(data.wide$Status)
# 把数据框内的变量放入矩阵内（Lasso回归法认不出数据框）
tmp.x <- model.matrix(~.,data.wide[c("Age","Sex_num","Status_num","TestA","TestB","TimeInterval_AC","TimeInterval_BC")])
tmp.y <- data.wide$TestC
# 进行Lasso回归
model.lasso <- glmnet(tmp.x, tmp.y, family="gaussian", nlambda=50, alpha=1, standardize=TRUE)
model.lasso #展示出一堆结果
# 寻找并展示lambda在第20行时，模型的信息
coef(model.lasso , s=model.lasso$lambda[20])
# 寻找并展示lambda等于0.076970时，模型的信息
coef(model.lasso , s=0.076970)
# 通过交叉验证，寻找最优模型
cv.model <- cv.glmnet(tmp.x, tmp.y, family="gaussian", nlambda=50, alpha=1, standardize=TRUE)
#得出交叉验证的模型图，竖线位置对应的红点的误差最小，对应下方的x轴位置为最优的lambda
plot(cv.model)
#得出最优lambda的值
cv.model$lambda.min
#寻找并展示lambda等于最优lambda的值时，模型的信息
coef(cv.model, s=cv.model$lambda.min)
## 3. 变量相对重要程度
#β值：自变量每变化一单位对因变量影响的大小。
#但是自变量本身尺度的不同，无法直接用β值表示自变量的重要性。
#所以提出标准化后的β，标准化后的β不受自变量本身尺度的影响；并且标准化系数绝对值的大小，可以作为评价自变量相对重要程度的参考指标
# 3.1 标准化回归系数：用于评估变量的重要性
install.packages('lm.beta')
library(lm.beta)
summary(model.full)#查看未标准化的β值
lm.beta(model.full)#查看标准化的β值
# 3.2 根据解释R平方的比例：对因变量y的方差解释的越多的自变量越重要
calc.relimp(model.full,type="lmg", rela=TRUE)#看Img法的结果（运行结果Img下的内容）
#得出通过交叉验证的图，这个图展示了各种变量的相对重要程度及相对重要程度的95%置信区间
boot <- boot.relimp(model.full, b = 1000, type = "lmg", rank = TRUE, diff = TRUE, rela = TRUE)
plot(booteval.relimp(boot,sort=TRUE))
# 3.3 根据Lasso保留变量的顺序
#下方x轴：log( λ）
#上方x轴：模型中变量个数
#先被排除的变量，相对不重要
plot(model.lasso, xvar="lambda", label=TRUE)
#上面的图如果看不清，可以逐个模型提取出来看，就知道被排除的是哪一个变量
# 展示只有1个变量是，回归模型长什么样
coef(model.lasso , s=model.lasso$lambda[2])
# 展示只有2个变量是，回归模型长什么样
coef(model.lasso , s=model.lasso$lambda[3])
# 展示只有3个变量是，回归模型长什么样
coef(model.lasso , s=model.lasso$lambda[5])
# 展示只有4个变量是，回归模型长什么样
coef(model.lasso , s=model.lasso$lambda[7])
# 展示只有5个变量是，回归模型长什么样
coef(model.lasso , s=model.lasso$lambda[12])
# 展示只有6个变量是，回归模型长什么样
coef(model.lasso , s=model.lasso$lambda[21])
# 展示只有7个变量是，回归模型长什么样
coef(model.lasso , s=model.lasso$lambda[34])
###第三章-7.多项式回归、分段回归、限制性立方样条、Lowess平滑方法
#操作前准备
rm(list = ls())
## install packages
install.packages('segmented')
install.packages('splines')
install.packages('ggplot2')
install.packages('Hmisc')
install.packages('rms')
library(segmented)
library(splines)
library(ggplot2)
library(Hmisc)
library(rms)
## load data
load('diagnosis.RData')
## add labels
data.wide$Status <- factor(data.wide$Status,
levels=c(0,1),
labels=c("No Disease", "Disease"))
data.wide$Sex <- factor(data.wide$Sex,
levels=c(0,1),
labels=c("Male", "Female"))
# attach data.wide
attach(data.wide)
# 0. 创建模型：model1 TestC~Age
model1 <- lm(TestC ~ Age)
y_model1 <- predict(model1)
# 1. 多项式回归（Polynomial Regression）
# 一元n次多项式：①V次:y=a+bx;②2次:y=a+bx+cx2;③3次:y=a+bx+cx2+dx3:④……
# 一般加入多项式最大为三次方，加入超过三次方会导致模型过度拟合的问题
# 生成多项式回归模型
model1.poly2 <- lm(TestC ~ Age + I(Age^2))
model1.poly3 <- lm(TestC ~ Age + I(Age^2) + I(Age^3))
# 查看生成的模型
summary(model1.poly2)
summary(model1.poly3)
# 比较两个模型
anova(model1,model1.poly2)#p＜0.05，说明两个模型有差异，后一个模型比前一个模型好
anova(model1.poly2,model1.poly3)#p＞0.05，说明两个模型无差异，不用加入三项式
# 预测y值
y_model1.poly2 <- predict(model1.poly2)
y_model1.poly2
y_model1.poly3 <- predict(model1.poly3)
y_model1.poly2
# 2. 分段回归（Segmented regression）
# 优点：分段回归是按照拐点分开的几段直线，不会像多项式那样牵一发而动全身
# 缺点：但是拐点处会比较生硬
# 某些模型，拐点处有一定的意义，但是我们操作中大多数时候更加希望得到光滑的曲线
model1.segmented1 <- segmented(model1)
summary(model1.segmented1)#Est下是拐点；Estimate下是截距和斜率
intercept(model1.segmented1)#查看每段直线的截距
slope(model1.segmented1)#查看每段直线的斜率#slope2的结果等于summary运行出两个斜率的和
#Age是自变量，TestC是因变量。pch=1指定了点的形状，cex=1.5调整了点的大小
plot(Age,TestC, pch=1, cex=1.5)
#添加一条直线，a是截距，b是斜率
abline(a=coef(model1)[1],b=coef(model1)[2],col="red",lwd= 2.5)
#添加分段线性模型的图形，col='blue'指定了线的颜色，lwd=2.5指定了线的宽度
plot(model1.segmented1, col='blue', lwd= 2.5 ,add=T)
#如果需要，可用psi指定拐点(35,50)后,再拟合分段线性模型
model1.segmented2 <- segmented(model1,psi=c(35,50))
#在图上添加另一个分段线性模型的图形，这个模型的断点是预先指定的（35和50）
plot(model1.segmented2, col='green', lwd= 2.5 ,add=T)
# 3. Spline样条回归
# Spline样条回归是结合多项式回归和分段回归的方法，可以灵活的展示自变量和因变量之间的关系
# 限制性立方样条是Spline样条回归最常用的一种，所以我们只学限制性立方样条
# 限制性立方样条两段是限制为线性的，这样可以，克服两端数据不足导致的过度拟合
# 在每个节点两侧，y的预测值相等(连续)，一阶导数相等(斜率相等)，二阶导数相等(平滑)
# 判断取几个节点更好
# Spline样条回归对节点数量不敏感，不需要担心节点数量选择会影响结果
model.spline1 <- lm(TestC ~ rcs(Age, 3))#选择3个节点
model.spline2 <- lm(TestC ~ rcs(Age, 5))#选择5个节点
AIC(model.spline1,model.spline2)#3个节点的模型AIC值比5个节点的更小，所以3各节点比5个节点更好
# 预测y值
y_model.spline1 <- predict(model.spline1)
y_model.spline2 <- predict(model.spline2)
# 另一种方法，使用ggplot画图
# 3个节点
ggplot(data=data.wide, aes(x=Age, y=TestC)) +
geom_point()+geom_smooth(aes(x=Age, y=TestC), method="lm", formula=y ~ rcs(x, 3), se=FALSE, colour="blue")
# 5个节点
ggplot(data=data.wide, aes(x=Age, y=TestC)) +
geom_point()+geom_smooth(aes(x=Age, y=TestC), method="lm", formula=y ~ rcs(x, 5), se=FALSE, colour="red")
# 把FALSE改成TRUE，可以显示95%置信区间
ggplot(data=data.wide, aes(x=Age, y=TestC)) +
geom_point()+geom_smooth(aes(x=Age, y=TestC), method="lm", formula=y ~ rcs(x, 5), se=TRUE, colour="red")
# 4. Lowess平滑方法（局部加权回归散点平滑法）
#原理：把自变量x从小排到大，以每一个数据点为中心，前后各截取一段数据作为加权的线性回归模型预测出对应于刚刚选取的x值的y的预测值，然后对每个x值重复这一段过程，得到所有x对应的y的预测值，连成一条直线，为Lowess平滑曲线
#优点：不仅可用于连续变量，还可用于分类变量
#优点：非常灵活，不需要对模型结构做假设，也不需要设置为二次项还是三次项方程
#缺点：没法得出具体的回归方程、决定系数等，只能给出一个图型，所以无法通过方程总结因变量与自变量的关系，没有外推性
#用处：通常用于在非常杂乱的数据探索因变量与自变量的关系，并通过肉眼可辨识的方式展示出来
#想进行下一步分析，可以将自变量分成多个组，利用多个组的y值的均值（Lowess回归的均值）作为该组的赋值，再利用这个赋值做回归
#用于连续变量
plot(Age,TestC, pch=1, cex=1.5)
lines(lowess(Age, TestC), col=2,lwd= 2.5)
#用于分类变量
plot(Age, Status, pch=1, cex=1.5)
lines(lowess(Age, Status), col=2,lwd= 2.5)
###第三章-8.广义线性模型（Generalized Liner Model (GLM)）
#一般线性模型有一个非常重要的假设：残差项必须服从正态分布
#但是当因变量y分布比较特殊的时候，残差正态分布很难被满足，这个时候需要使用广义线性模型
#广义线性模型的右侧：为自变量的线性组合，与一般线性模型一样
#广义线性模型的左侧：连接函数g，使得g(E[Y])等于右侧线性组合
## 广义线性模型要求因变量Y服从指数分布族：
# 连续变量 —— 正态分布 —— 一般线性回归
# 二分类变量 —— 二项分布 —— Logistic回归
# 计数数据 —— 泊松分布 —— 泊松回归
# 操作前准备
rm(list = ls())
## install packages
install.packages('ggplot2')
install.packages('Hmisc')
install.packages('MASS')
install.packages('rms')
install.packages('pROC')
library(ggplot2)
library(Hmisc)
library(MASS)
library(rms)
library(pROC)
## load data
load('diagnosis.RData')
## add labels
data.wide$Status <- factor(data.wide$Status,
levels=c(0,1),
labels=c("No Disease", "Disease"))
data.wide$Sex <- factor(data.wide$Sex,
levels=c(0,1),
labels=c("Male", "Female"))
# new variables
data.wide$TestA_log <- log(data.wide$TestA)
data.wide$TestB_log <- log(data.wide$TestB)
data.wide$TestC_log <- log(data.wide$TestC)
# attach data.wide
attach(data.wide)
# 1. Logistic回归
# 原理：求对于给定的自变量X，因变量Y=1的概率（p（X）=Pr(Y=1丨X)）
# 用logit函数作为连接函数(log(p(X)/(1-p(X)))=β0+β1X)
# 建立Logistic回归模型
# family=binomial()是告诉R软件要用哪一种指数分布族，下面的是二项分布
logistic1 <- glm(Status ~ Age,
data=data.wide,
family=binomial())
# 查看结果
summary(logistic1)#看截距和斜率（β值）
confint.default(logistic1)#看截距和斜率（β值）的95%置信区间
# 求OR值
exp(coef(logistic1))
# 预测概率值
p.logistic1 <- predict(logistic1, type = "response")
p.logistic1
# 求自变量（年龄）与因变量（概率预测值）的关系曲线
ggplot() + geom_line(aes(x = Age, y = p.logistic1))
# 做变量筛选
# 构建所有自变量（其中包括log函数计算的变量）与因变量的模型
logistic.full <- glm(Status ~ Age + relevel(Sex,"Male") + log(TestA) + log(TestB) + log(TestC), data=data.wide, family=binomial())
# 用AIC法筛选变量
logistic.select <- stepAIC(logistic.full)#最后筛选出两个变量：log(TestB)、log(TestC)
# 查看筛选结果
summary(logistic.select)
# 用新的模型预测概率值
p.select <- predict(logistic.select, type = "response")
# 2. ROC曲线（ROC curve）
# 进行ROC分析
roc.logistic1 <- roc(Status,p.logistic1)#后面出现的Setting levels、Setting direction并不是运行失败，而是pROC包设定Status变量的两个水平被为“control”（对照组，即“无疾病”）和“case”（病例组，即“有疾病”），这是pROC包在进行ROC分析时，对二分类结果变量的水平进行的命名。并且pROC包自动设定了比较的方向为controls < cases，即对照组的预测概率小于病例组的预测概率。
roc.logistic1#“Area under the curve:”后是ROC曲线下面积的值
# 看ROC曲线下面积的95%置信区间
ci(roc.logistic1)
# 绘制ROC曲线
plot(roc.logistic1)
# 比较ROC曲线下面积(AUC)
# 在做另一个ROC分析
roc.select <- roc(Status,p.select)
# 把另一个ROC分析曲线加到图上
plot(roc.select, add=TRUE,col=2)
# 比较两个ROC曲线下面积是否有不同
# 运行后可以看到两个面积差距较大，但是p值没差异
# 所以我们不能仅仅通过p值判断两个曲线下面积是否有差异
roc.test(roc.logistic1 , roc.select, method="delong")
# 确定最优阈值
plot(roc.select)
# 求阈值
# “best.method="youden”准则：让敏感度和特异度的和最大
coords(roc.select,"b",best.method="youden")# threshold下显示最佳阈值
# “closest.topleft”准则：选离左上角最近的点做阈值
coords(roc.select,"b",best.method="closest.topleft")
# 3. 列线图（Nomogram）
# 列线图是风险得分计算器+风险得分与概率对照表
# 主要目的是展示预测模型和作为计算工具
# 如果模型的目的不是预测，不建议使用列线图
#这行代码创建了一个datadist对象，它包含了data.wide数据集的分布信息
#这个对象会被用来设置rms包中模型的默认分布
ddist <- datadist(data.wide)
#这行代码设置了全局选项datadist，使得rms包中的模型都会使用之前创建的ddist对象
options(datadist='ddist')
#建立模型
logistic.lrm <- lrm(Status ~ Age + Sex + TestA_log + TestB_log + TestC_log, data=data.wide)
#创建列线图
#fun=plogis表示使用logistic函数（sigmoid函数）来转换线性预测值
#lp=F表示不使用对数尺度（log-scale）
#funlabel="Disease Risk"为转换函数设置了标签
nom.full <- nomogram(logistic.lrm, fun=plogis,lp=F, funlabel="Disease Risk")
#画图
plot(nom.full)
#列线图的阅读方法：把每一项变量对应的值对应上面的分数，把每一个变量得出的分数相加，找到总分对应下面的患病风险的值
扫一扫
18万+
11-23
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
