第一课：安装与基本操作

R的扩展包在R官网CRAN；另外，R官网还包含很多扩展资料，包括源代码，手册，FAQ，推荐书籍等，该网站本身就是学习R非常好的资源站点；缺点：不够规范，不容易上手，需要付出较大努力；其次，R的扩展包太多了，需要很多时间查找和学习，很多扩展包的学习比学习R本身还要困难；

1 很多时候需要对原始文件进行格式转换才可以用，比如去除缺失值，调整行和列的内容等；首先使用table函数进行样品频数统计；aggregate函数可以进行频数统计；

R的32位和64位版本是同一个安装包；

最新版的Win10内置了Linux shell，可以在里面安装R；MAC的OS系统自带Unix终端，也可以安装R；linux主要可以使用R源代码来安装R；一般版本，奇数版本为开发版，偶数版本为稳定版；

2 R软件的运行与设置

Win的32位和64位是同一个安装包，虽然32位版本快一些，但是只能支持3G的内存，64位对内存没有限制；列出所有选项可以用ls命令来完成；程序包可以对R的扩展包进行操作；帮助菜单有很多获取帮助的内容，使用R软件时，获取帮助也非常重要；

3 Mac安装R后，直接在终端运行，大写R 即可进入R ; 完全面向对象，敲getwd命令即getwd（）会返回当前工作目录；运行max命令挑选最大值；一定要将括号补齐才是完整的表达式；如需反复使用，可以将代码储存在.R的文件中反复调用，这就是R的脚本；运行R脚本可以使用R script命令，或者将运行的R保存为R工作空间印象，退出时选择y保存工作印象会在退出时的工作目录中生成一个Rdata的文件，双击打开该文件即可运行R

4 通过命令行设置R

启动R时，R会加载Rprofile.set文件，可通过修改它来设置界面R；

5 R studio作为R的集成开发环境，本身不包含R软件；

R studio安装时会自动搜索R，直接打开即可使用；

左上为代码窗口用来打开脚本文件；左下为R控制台；右上为环境和历史记录，右下为功能窗口，包括环境、功能、帮助文档窗口等；

File-new file-script 则是脚本文件

选中几行代码，直接run则是运行当前行；run旁边的rerun是重复运行，再右侧运行脚本中所有代码；source with echo会显示具体的脚本信息，包括脚本运行的选项参数；荧光棒检查代码，放大镜有查找和替换功能，保存图表可将代码保存为不同文件类型，打开不同文件显示也会有差别；通常在控制台中直接敲代码，一万行入门！

6 基本操作

R中函数命令必须加括号，因为常量、变量、函数、图形，在R中都是对象，都可以用自定义的字母来表示；

Getwd函数获取当前工作目录，tools-global options可以修改选择长期工作路径,修改后list.files命令可以查看目录下包含的文件,dir命令可以完成同样的工作；R中直接以字符定义变量，但是变量名不可以数字开头,符号之间不可有空格；也可用等号进行赋值,但进行假设检验时会与其中等号混淆，所以不推荐用等号赋值；option加连字符是赋值快捷键；连字符和小于号是向右赋值，反过来即可；<<-表示强制赋值给全局变量而不是局部变量；

Mean函数对第一个参数求均值，z<-mean(c(1,2,3,4,5))即可得到这5个数的均值;

查看当前工作空间变量及函数：

Ls（）得到当前变量/对象；ls.str(得到具体信息)；单独str函数则可以得到每个变量的详细信息；ls(all.names=TRUE)

删除对象或值：rm(),删除所有对象，首先将所有对象赋给list，然后直接rm（list）即可删除全部对象；

上下光标选择近期命令；history函数列出历史记录，history（25）可列出最近25条历史记录；

命令过多-清空屏幕：control+L组合键

为防止死机定期保存工作空间：save.image(),但只能保存当前数据和函数，图片需另外保存，q退出软件

7 R包安装

R-cran-task view浏览各类R包；

专门处理生物信息的包：bioconductor

install.packages(“包”)自动下载安装；

已安装vcd包；

.libpaths()可以显示库所在的位置；

绝大部分可以线上安装并且自动安装所需的依赖包；

源代码安装示例：安装ABCanalysis包，该包依赖plotrix包；

update.packages（）可以更新安装的软件包

R包使用

library()展示库中的包；

library（vcd）载入vcd包，此时不用加引号；require（）也可加载包

R自身的包 提供了一些默认函数和数据集；

函数search（）可得知哪些包已经加载并可使用；

使用包vcd举例，help（package=“vcd”）可得到该包的帮助文档，包含描述文件，用户指导手册，代码展示，实例展示,library(help=“vcd”)可以查看包的信息，包含一些包的基础内容；index下面就是数据集，一般是给包举例使用的;

ls("package:vcd”)可以显示vcd包中包含的所有的函数

data(package="vcd”)可以显示该包中所有的数据集

detach("package:vcd") 移除该包

remove.packages(“vcd”)可以从硬盘上彻底删除该包

R包批量移植：

installed.packages()

installed.packages()[,1]得到包的名字，方法：使用下标访问数据框的第一列；

Rpack <- installed.packages()[,1]

save(Rpack,file="Rpack.Rdata”)将R包的名字保存成一个文件,将此文件copy到另一台设备，可以用load函数来加载；

for(i in Rpack) install.packsges(i) 将旧机器所有的R包安装在新设备上；如果重复没关系，R会跳过已安装的扩展包；

8 获取帮助

help.start()可以打开并查阅相关的帮助文档；

help("sum”)也可以显示sum函数的帮助文档；

？+函数名也可以出现帮助文档；

Args（函数名可以出现函数的主要功能和参数）

Example（函数名）可以出现函数示例，包括普通函数和绘图程序，都有示例；

demo(graphics)会列出用R绘制的图

help(package=ggplot2) 查询安装的包的名字

Vignette（“”）为另一种形式的文档，包括简介，教程，开发文档等；

已安装但使用help函数没有办法调取帮助文档，此时需要先用library函数载入包再调取帮助文档；

??ggplot2  不加载包也可以出现帮助文档；

help.search

help.search("heatmap”)可以搜索到自己安装的或者stats中的heatmap函数用于绘制热图；与??heatmap的功能等同；

apropos("sum”)可以列出所有包含sum关键字的内容

apropos("sum",mod="function") 找出包含sum的内容并只显示函数；

RSiteSearch("leukemia") 该函数可使用默认浏览器访问R官网，搜索文档；

google上Rseek.org网站内是只与R有关的内容；

9  Excel 案例

command+1 快捷键打开单元格格式

R中也需准确分辨数据类型；

R中一般行是观测值，列是变量值，每列的数据类型一般一致，为了方便计算；

###字符串

nchar("hello world") # nchar函数统计字符串的长度，包括空格

month.name

nchar("month.name") #记住，字符串一定要加引号

nchar(month.name) #记住，字符串一定要加引号

length(month.name) #length函数给出的是元素个数，nchar函数给出的是字符串的长度也就是有几个字母和空格

nchar(c(12,2,345))  #此时nchar返回的也是字符串长度，这是是将数字向量转化为字符串向量来处理

paste(c("everybody","loves","stats"))  #paste函数可将多个字符串合并为一个

paste("everybody","loves","stats",sep="-")  #paste函数可将多个字符串合并为一个，这里不要用c，否则出来的字符串还是有双引号的，默认使用空格分隔，可以使用sep函数设定分隔符，

names <- c("moe","larry","curly")

names

paste(names,"loves stats") #向量和字符串相加，并不是简单的加到之后，而是对names中每个元素分别处理

#输出的是"moe loves stats"  "larry loves stats" "curly loves stats"

substr(x=month.name,start=1,stop=3) #substr为提取字符串，该代码可以提取出每个月份名字前3个字母

#结果为[1] "Jan" "Feb" "Mar" "Apr" "May" "Jun" "Jul" "Aug" "Sep" "Oct" "Nov" "Dec"

temp <- substr(x=month.name,start=1,stop=3)

toupper(temp)  #toupper函数可以将字符串转化为大写

tolower(temp)  #将此转化为小写

#如果在此基础上将首字母变为大写，可以使用sub或者gsub函数，sub是单个替换，gsub是全局替换，此处需要用到正则表达式

gsub("^(\\w)","\\U\\1",tolower(temp)) #倒三角表示首字母，短斜杠w表示字符集简写，代表所有小写字符，短斜杠U表示转化为大写，1代表只转换一次，tolower（temp）代表转换的对象

#上面结果是不行的，只给每个前面加了一个U，所以要设置参数，perl=1，支持perl类型的正则表达式

gsub("^(\\w)","\\U\\1",tolower(temp)，perl=T) #倒三角表示首字母，短斜杠w表示字符集简写，代表所有小写字符，短斜杠U表示转化为大写，1代表只转换一次，tolower（temp）代表转换的对象,perl=T代表支持perl类型的正则表达式

gsub("^(\\w)","\\U\\1",tolower(temp),perl=T) #倒三角表示首字母，短斜杠w表示字符集简写，代表所有小写字符，短斜杠U表示转化为大写，1代表只转换一次，tolower（temp）代表转换的对象,perl=T代表支持perl类型的正则表达式

gsub("^(\\w)","\\L\\1",tolower(temp),perl=T) #切换首字母小写，将U换成L即可

?grep

#grep函数搜索字符串

x <- c("b","A+","AC")

x

grep("A+",x,fixed=T) #fixed参数指匹配A+的字符串所在的位置

grep("A+",x,fixed=F) #fixed参数指匹配A+的字符串所在的位置，如果此时是F，则根据正则表达式，A+指的是A到正无穷大，则AC也满足条件，会出现A+和AC两个字符串的下标，即2和3

match("A",x) #match函数可以进行字符串匹配，也会给出相应下标，但是不支持正则表达式，没有grep函数强大

match("A+",x)

#strsplit函数进行字符串的分割

path <- "/user/local/bin/R"  #是linux下的一个路径，用斜线/分割，strsplit函数则需要输入需要分割的字符串和分隔符作为参数进行分割

strsplit(path,"/")  #注意，此处的分隔符需要加双引号

#strsplit返回的值是列表而不是向量，这样是为了方便处理多个字符串，如果我们同时处理两个，则返回向量就不合适了，此处返回的是列表

strsplit(c(path,path),"/")  #注意，此处的分隔符需要加双引号

#fixed参数是为了控制是否使用正则表达式

#生成字符串成对组合的技巧，比如两个字符串生成他们所有的组合，也叫做笛卡尔集，需使用outer函数

face <- 1:23

suit <- c("spades","clubs","hearts","diamonds")

face

suit

outer(suit,face,FUN=paste)  # 内包含两个向量，参数FUN接一个函数，因为要连接字符串，所以使用paste函数

#更正，上面face应该是1：13，这样可以生成一副扑克牌，以下重复一遍

face <- 1:13

outer(suit,face,FUN=paste,sep="-") #使用sep参数指定连接符号

savehistory("~/Desktop/R-data/R-字符串处理.Rhistory")

2021-06-17

data3

#cbind和rbind也可用于matrix，但是两者使用时必须要求具有相同的行数或者列数，否则就会出现问题，测试一下

head(cbind(USArrests,state.division) ,20)

data3 <- head(cbind(USArrests,state.division) ,20)

data2

rbind(data3,data2)

#因为data3和data2的列数不一样，所以无法用rbind合并到一起

data1 <- head(USArrests,30)

data2 <- head(USArrests,10)

rbind(data1,data2)  #data1和data2有10行是重复的，将两者合并后会怎样

#在列重复的情况下，合并后相同的列会自动合并成一列，相同的行则会重复出现。这是完整的合并，没有去除重复项

data1 <- head(USArrests,30)

data2 <- tail(USArrests,30)

data3 <- rbind(data1,data2)

data3

rownames(data3)

length(data3)

length(rownames(data3))

##这里可见如果行重复了，合并都会保留，不会去掉重复项，但这样有时不便分析；R中操作相当于取子集，先获得没有重复项的行名，然后根据ID索引数据即可

duplicate(roenames(data3)) #先用duplicate函数判断是否有重复值

duplicated(roenames(data3)) #先用duplicate函数判断是否有重复值

duplicated(rownames(data3)) #先用duplicate函数判断是否有重复值

duplicated(data3) #先用duplicate函数判断是否有重复值

data4 <- data3[duplicated(data3),] #这样取出了重复部分，可以看到R已经将重复ID添加了1

data4

data5 <- data3[!duplicated(data3),] #这样取出了重复部分，可以看到R已经将重复ID添加了1,以上是取出重复部分，如果在索引中加上感叹号取反，就是取出非重复部分

data5

length(rownames(data3[!duplicated(data3)),]) #统计非重复部分有多少行

length(rownames(data3[!duplicated(data3)),]) #统计非重复部分有多少行

length(rownames(data3[!duplicated(data3)],))

rowanems(data5)

rownames(data5)

length(rownames(data5))

length(rownames(data3))

length(rownames(data3)) #已知data3有重复的部分，使用unique函数可以一步去重

data6 <- unique(data3)

length(rownames(data6))

length(rownames(data6)) #使用unique函数去重之后得到data6，然后计算data6的行名长度是50，可知data6是已经去重的数据

###数据转换第三课

sractm <- t(mtcars #t函数实现数据框反转，行列相互变换

）

sractm <- t(mtcars） #t函数实现数据框反转，行列相互变换

?t

sractm <- t(mtcars)

sractm

letters#反转单独一行用rev函数，rev函数针对的对象是向量

rev(letters)

rev函数反转数据框的操作，以下是演示

#rev函数反转数据框的操作，以下是演示

women

rownames(women)  #使用rownames函数访问行名

rev(rownames(women) ) #使用rownames函数访问行名

rev(rownames(women) ) #rev

rev(rownames(women) ) #rev函数对行名进行反转

women[rev(rownames(women) ) ,]#用翻转后的行名作为索引

women[rev(rownames(women) ) ,]#用翻转后的行名作为索引，就可以得到完全反过来的women数据框

data7 <- women[,height]#修改数据框中的值，将height单位英寸转化为厘米，需要将height中每个值*2.54

data7 <- women[,"height"]#修改数据框中的值，将height单位英寸转化为厘米，需要将height中每个值*2.54

data7

data8 <- data7*2.54

data8

women1 <- data.frame(height,women[,"weight"])

women1 <- data.frame(data8,women[,"weight"])

women1

#上面的操作出来的数据不太对

women1 <- data.frame(height=data8,weight=women[,"weight"])

women1

##上面两个代码可以运行，但是都是计算后需要修改，这样效率低，如果需要修改的值很多，那么更不高效

data9 <- transform(women,height=height*2.54) ##为避免麻烦可以使用transform函数，该函数可以非常方便修改数据框中的值,一步到位

data9

data9 <- transform(women,cm = height*2.54)#也可以不在原数据上修改，多一列，就多一个名字而已

data9

sort(rivers)#数据框的排序，R中排序的函数主要有3个，sort，order和rank；sort是对向量进行排序

sort(state.name) #默认数字是按数字大小的顺序，字符串是按首字母字母表顺序

rev(sort(rivers)) #配合rev是对排序后的结果进行反转

mtcars[sort(rownames(mtcars))],]

mtcars[sort(rownames(mtcars)),]]

mtcars[sort(rownames(mtcars)),]

rownames(mtcars)#虽然sort不能对数据框排序，但是可以通过去取出单个列，排序之后，对排序结果进行索引的方式对整个数据框进行排序

sort(rownames(mtcars))

mtcars[sort(rownames(mtcars)),]

##以上是一步步写的步骤

sort(rivers#sort返回的是排序后的结果，但是order函数返回的是排序后的位置，举例

sort(rivers)#sort返回的是排序后的结果，但是order函数返回的是排序后的位置，举例

sort(rivers)

order(rivers)

river[9]

rivers[9]

#order函数的好处是因为返回的是索引值，所以可以直接用索引值访问数据，这样间接地对数据进行了排序

order(mtcars$mpg) #用order函数对mpg这一列的大小进行排序

mtcars$mpg

mtcars

mtcars$disp

order(mtcars$disp)

mtcars[order(mtcars$disp),]

sort(mtcars$disp)

mtcars[sort(mtcars$disp),]

mtcars[,sort(mtcars$disp)]

#据我的观察，虽然order提供的是索引值，但是对整个数据框排序时还是按照从小到大的顺序进行排序

#order能对非首列进行排序，但是sort只能对第一列进行排序

mtcars[order(-mtcars$disp),]#如果要使用逆向排序，只需要在mtcars前面加上减号即可

#rank和其他两个函数不一样，是求秩的函数，返回的是该向量对应元素的排名

##数据转换4

#对数据框进行数学计算

WorldPhones #全球8个地区7年中的手机数量

WP <- as.data.frame(WorldPhones)#求每一年全球总电话呼叫数量和7年每个大洲平均呼叫数量,首先将矩阵转化为数据框

WP

rs <- rowSums(WP)

rs

cm <- colMeans(WP) #colmeans计算每个大洲的平均数

cm

#其实rowsums和colmeans就是对行取和，对列取平均

WP1 <- data.frame(WP,rs)

WP1

WP2 <- data.frame(cm,WP1)

WP2

#添加一列比较方便，添加行这里出了问题

WP2 <- cbind(WP,rs)

WP2

WP2 <- rbind(WP2,cm)

WP2

rowname <- rownames(WPA)

rowname <- rownames(WP2)

rowname

rowname[8] <- mean()

rowname[8] <- "mean"

rowname

rownames(WP2) <- rowname

WP2

#以上的操作也将mean和rs值添加进去了

WP

sum <- rowSums(WP)

mean <- colMeans(WP)

sum

mean

WP2 <- cbind(WP,sum)

WP2

WP3 <- rbind(WP2,mean)

WP3

rownames(WP3)[8] <- mean

rownames(WP3)[8] <- "mean"

WP3

WP3 <- cbind(WP,sum)

WP3

WP3 <- rbind(WP3,mean=mean)

WP3

#但是注意，右下角这个值并不是sum的总和，这个值被重复使用了

sum1 <- WP3[,sum1]

sum1 <- WP3[,sum]

sum1 <- WP3[,"sum"]

sum1

sum1[66747.57] <- NA

sum1

sum1[66747.57] <- "NA"

sum1

sum1$66747.57 <- NA

sum1

?sum1

??sum1

str(sum1)

sum1 <- numberic(sum1)

?numberic

sum1 <- as.numeric(sum1)

sum1 <- sum1[-1]

sum1

#为了避免麻烦，可以使用apply系列函数

#apply：操作对象可以是矩阵，数组或者数据框，margin是维度的下标，分为1和2，1代表对行进行处理，2代表对列进行处理，FUN代表函数，可以使用sum，mean等进行计算

WP

#以上的WP是对WorldPhones的简称

apply(WP, MARGIN = 1, FUN = sum) #该代码表示对WP数据框中的行进行求和

apply(WP, MARGIN = 2, FUN = mean) #该代码表示对WP数据框中的每列求平均值

apply(WP, MARGIN = 2, FUN = log) #FUN参数也可以计算方差或者log值

?lapply

#l代表list，返回的是列表，s意思是simplify，简化后的，返回向量或矩阵，比列表要简单

state.center #该数据集是列表，不可以使用普通的apply函数，可以使用lapply函数

lapply(state.center,FUN=length) #列表没有行列之分，所以省去margin这个词参数

lapply(state.center,FUN = sum)

lapply(state.center,FUN = mean)

lapply(state.center,FUN = log())

lapply(state.center,FUN = log

)

class(lapply(state.center,FUN = mean))

sapply(state.center,FUN = mean)

sapply(state.center,FUN = sum)

class(sapply(state.center,FUN = sum)) #sapply函数同样可以操作，但这次返回的是向量值

class(sapply(state.center,FUN = length) #sapply函数同样可以操作，但这次返回的是向量值

class(sapply(state.center,FUN = length) ）#sapply函数同样可以操作，但这次返回的是向量值

class(sapply(state.center,FUN = length))

sapply(state.center,FUN = length)

sapply(state.center,FUN = sum)

#tapply稍有不同，它处理因子数据，根据因子分组，然后对每组分别处理。第一个参数是数据集，第三个参数是要使用的函数，最重要的是第二个参数也就是index，必须是因子类型，利用该因子对第一个参数内的数据进行分组

state.name #美国50个洲名字

state.region #对50个洲进行分区

state.division #division才是对50个洲进行分区，分为了几个大区

#统计各大区包含哪几个州

tapply(state.name,state.division,FUN = length) #第一个参数是要处理的数据集，第二个参数是因子，也就是给name归类找了一个参照，第三个是函数，也就是要看每个因子内元素的长度

tapply(state.name,state.region,FUN = length)

#apply系列中最为强大的是FUN参数，包含很多函数，可以进行很多种计算

#R中的数据中心化和标准化处理 中心化：数据集中各数据减去数据均值；标准化：中心化后除以数据集的标准差；目的：消除量纲对数据结构的影响

heatmap(state.x77) #由于state.x77各数据差别太大，所以绘制热图之后没有意义

x <- c(1,2,3,6,3) #举例说明

mean(x)

x-mean

x-mean(x)

sd(x)

y <- sd(x)

z <- x-mean(x)

bz <- x/y

bz

bz <- z/y

bz

scale(x,center = T,scale = F)#也可以直接使用scale函数进行中心化和标准化

scale(state.x77,scale = TRUE)

y <- scale(state.x77,scale = TRUE)

heatmap(y)

q()

reshape2

##2021-06016 reshape2

#reshape2包对数据进行格式转换

x <- data.frame(k1)

x <- data.frame(k1=c(NA,NA,3,4,5),k2=c(1,NA,NA,4,5),data=1:5)

x <- data.frame(k1=c(NA,2,NA,4,5),k2=c(NA,NA,3,4,5),data=1:5) #设置x和y两个数据集

x

y

#x和y两个数据集的内容相互补，但是由于是两个数据集，所以不能用cbind和rbind函数，否则会连成一片，分不清数据来源，用merge函数对两数据集中的公有变量进行合并，下面分别依据看

merge(x,y,by="k1")#分别依据k1,k2进行合并

y <- data.frame(k1=c(NA,2,NA,4,5),k2=c(NA,NA,3,4,5),data=1:5) #设置x和y两个数据集

y

merge(k1,k2,by)

merge(k1,k2,by="k1")

merge(x,y,by="k1")

merge(x,y,by="k2")

merge(x,y,by="k2",incomparable=T) #incomparable选项表示丢掉NA的情况

merge(x,y,by = c("k1","k2"))

#但是merge函数比较麻烦，使用reshape包可以将数据转换成任何需要的形式

#但是merge函数比较麻烦，使用reshape包可以将数据转换成任何需要的形式?

install.packages("reshape2")

library(reshape2)

airquality

airquality #内置数据集，纽约1973年5-9月每天空气质量情况

head(airquality)

names(aiqqu)

names(airquality) <- tolower(names(airquality)) #使用tolower函数将列名改成小写

head(airquality)

melt(airquality)

#对该数据集融合之后每一行都是独立的标志符和变量的组合，没有重复项，融合之后数据变成三列，分别是序号，指标，指标数值，variable是一个因子类型，这是一个把宽数据变成长数据的操作

#但是需要告诉melt函数还需要month和day作为行的观测；也就是规定好哪部分作为行的观测，哪部分作为列的观测值，ID的参数是vars

aql <- melt(airquality,vars=c("month","day"))

aql <- melt(airquality,id.vars=c("month","day"))

#更正，id.vars才是规定行的观测的参数

head(aql)

#这下就正确了

#melt是融合，后面还需要重铸，用cast函数，包括dcast和acast，处理数据框用dcast函数，acast函数返回向量、矩阵、数组

aqw <- dcast(aql,month+day ~varibale) # 这里用dcast重铸数据，month+day~varibale表示以这两个作为主要的，波浪线表示观测值

aqw <- dcast(aql,month+day ~variable) # 这里用dcast重铸数据，month+day~varibale表示以这两个作为主要的，波浪线表示观测值

aqw

aqw <- dcast(aql,month ~variable) # 只显示每个月的情况

aqw <- dcast(aql,month ~variable,fun.aggregate = mean) # 只显示每个月的情况，此时需要指定对每天的数据如何处理，fun.aggregate

aqw

aqw <- dcast(aql,month ~variable,fun.aggregate = mean,na.rm=TRUE) # 只显示每个月的情况，此时需要指定对每天的数据如何处理，fun.aggregate,使用na.rm参数去除缺失值

AQW

aqw

#把reshape里的dcast和acast的案例都运行一遍就可以熟悉reshape2包的使用

？dcasst

?dcast

#Chick weight example

ChickWeight

names(ChickWeight) <- tolower(names(ChickWeight))

chick_m <- melt(ChickWeight,id=2:4,na.rm=TRUE)

chick_m

head(ChickWeight)

head(chick_m)

dcast(chick_m,time=variable,fun.aggregate = mean)

dcast(chick_m,time=variable, mean)

dcast(chick_m,time~variable, mean)

dcast(chick_m,dir.exists()~variable, mean)

dcast(chick_m,diet~variable, mean)

dcast(chick_m,diet~time, mean)

install.packages("tisyr")

install.packages("tidyr")

install.packages("dplyr")

library(tidyr)  #tidyr包用来处理tidy data，也就是整洁，干净的数据

help("tidyr")

help(package="tidyr")

#tidyr包有4个重要函数  ；gather将宽数据转换为长数据；spread将长数据转化为宽数据；unite将多列合并为一列；seaprate将一列拆分为多列

tdata <- mtcars(1:10,1:3)

tdata <- mtcars[1:10,1:3]

tdata

#以mtcars作为演示，这是一个明显的tidy data

tdata <- data.frame(tdata,name=rownames(tdata))

tdata

tdata <- data.frame(name=rownames(tdata),tdata)

tdata

tdata <- [,-"name.1"]

tdata <- tdata[,-"name.1"]

tdata <- mtcars[1:10,1:3]

tdata

tdata <- data.f

tdata <- data.frame(name=rownames(tdata),tdata)

tdata

gather(tdata,key="K",value="Value",cyl,disp,mpg)

gather(tdata,key="K",value="Value",cyl：disp,mpg) #添加冒号参数可以指定将哪些列合并到一起

gather(tdata,key="K",value="Value",cyl:disp,mpg) #添加冒号参数可以指定将哪些列合并到一起

gather(tdata,key="K",value="Value",cyl:disp,-mpg) #减号指不需要转换的列，这样可以将mpg单独放到一列中

gather(tdata,key="K",value="Value",cyl:mpg,-disp) #这样可以将disp单独放到一列中

tdata

gather(tdata,key="K",value="Value",2:4)#为避免敲列的名字出错，可以直接使用列的编号，2：4也就是将2-4列放在一列中

##注意：列名和列编号都可以使用

#相比于其他函数，gather可以让固定列不变，而其他列进行转换

tdata <- gather(tdata,key="K",value="Value",2:4)

tdata

gdata <- gather(tdata,key="K",value="Value",2:4)

gdata

gdata <- tdata

gdata

spread(gdata,key = "Key",value = "Value") #key代表观测，value代表观测值，这样直接将3个观测拆开了

spread(gdata,key = "K",value = "Value") #key代表观测，value代表观测值，这样直接将3个观测拆开了

sdata <- spread(gdata,key = "K",value = "Value") #

sdata

unite(sdata,key = "K",value = "Value") #unite函数与spread函数相对应，spread表示拆分，unite函数表示合并列

#以上一行代码不算，接下来介绍的是separate函数

df <- data.frame(c=c(NA,"a.b","a.d","b.c"))

df

separate(df,col=x,into = c(A,B))

separate(df,col=x,into = c("A","B"))

df <- data.frame(x=c(NA,"a.b","a.d","b.c"))

df

separate(df,col=x,into = c("A","B"))

separate(df,col=x,into = c("A","B")) #直接依据点进行了拆分处理

df <- data.frame(c=c(NA,"a-b","a-d","b-c"))  #刚刚没加的参数是分隔符，因为会默认识别分隔符，将df中的分隔符替换成连字符

df

？separate

?separate

separate(df,col=x,into = c("A","B"),sep = "-") #依据连字符进行了拆分处理

separate(df,col=c,into = c("A","B"),sep = "-") #依据连字符进行了拆分处理

df

sdf <- separate(df,col=c,into = c("A","B"),sep = "-")

sdf

unite(sdf,col = "AB",sep = "-") #unite对separate后的结果进行合并，sdf是需要操作的数据框，col是新的列，AB是列名，sep是新建的分隔符

##dplyr包

library(dplyr)

ls("package:dplyr") #dplyr包功能太过强大，使用ls查看dplyr包中的函数，dplyr包不仅可操作性单个表，还可对不同的表进行操作

help(package = "dplyr")

iris

?filter

dplyr::filter(iris,Sepal.Length >7) #使用filter函数，根据sepal.length>7这个条件进行过滤

dplyr::filter(iris,Sepal.Length <8)

df1 <- data.frame(iris[c(10,15),]）

df1 <- data.frame(iris[c(10,15),])

df1

df1 <- data.frame(iris[c(10,15),]) #取iris第10行和15行构建一个数据框

df1 <- data.frame(iris[c(1：10,1：15),]) #取iris的1-10行和1-15行构建一个数据框

df1 <- data.frame(iris[c(1:10,1:15),]) #取iris的1-10行和1-15行构建一个数据框

df1

disdf1 <- dplyr::distinct(df1)  #distinct函数去重，dplyr函数众多，可能与其他函数重名，所以使用双冒号dplyr::distinct函数直接使用该函数，不会出错

disdf1

dplyr::slcie(iris,10:15) #slice函数意为切片，可以取出数据的任意行，10：15指的是取出10-15行

dplyr::slice(iris,10:15) #slice函数意为切片，可以取出数据的任意行，10：15指的是取出10-15行

sdf <- dplyr::slice(iris,10:15)

sdf

rownames(sdf) <- c(10:15)

sdf

dplyr::filter(sdf,Sepal.Width>4)

dplyr::filter(sdf,Sepal.Width>3.5)

samdf <- dplyr::sample_n(iris,10) #sample_n函数随机取样，对iris随机抽取10行

samdf

samfracdf <- dplyr::sample_frac(iris,0.1) #sample_frac用于按比例随机选取

samfrac

samfracdf

#sample_frac的具体用法没讲

dplyr::arrange(iris,Sepal.Length) #arrange函数用于排序，第二个参数是排序依据，这里指按照花萼长度排序

dplyr::arrange(iris,desc(Sepal.Length) #arrange函数用于排序，第二个参数是排序依据，这里指按照花萼长度排序

dplyr::arrange(iris,desc(Sepal.Length)) #arrange函数用于排序，第二个参数是排序依据，这里指按照花萼长度排序,desc对相反方向进行排序

dplyr::arrange(iris,desc(Sepal.Length)) #arrange函数用于排序，第二个参数是排序依据，这里指按照花萼长度排序,desc对相反方向进行排序

?select

library(tidyverse)

iris <- as.tribble(iris)

iris <- as.tibble(iris)

iris <- as_tibble(iris)

iris

summarise(iris,avg=mean(Sepal.Length)) #还有统计学功能，这里就是对sepal.length这一列进行统计，计算花萼的平均长度

#该函数的好处在于可以直接取出一行进行计算

summarise(iris,avg=sum(Sepal.Length))  #这是计算总长度

summarise(iris,avg=first(Sepal.Length))

summarise(iris,avg=last(Sepal.Length))

iris

tail(iris)

#链式操作符%>%,两个百分号中间夹一个大于号，用于实现将一个函数的输出传递给下一个函数，作为下一个函数的输入

#在R studio中可以使用ctrl+shift+M快捷键输出

head(mtcars,20) %>% tail(mtcars,20)

head(mtcars,20) %>% tail(10)

#注意链式操作符是对前一处理的结果进行操作而不是对一个操作对象施加两个处理

dplyr::group_by(iris,Species) #group_by函数可以对数据进行分组

iris %>% group_by(Species) %>% summarise() #先用species进行分组，然后进行分组统计

iris %>% group_by(Species) %>% summarise(avg=mean(iris)) #先用species进行分组，然后进行分组统计

iris %>% group_by(Species) %>% summarise(avg=mean(Sepal.Length)) #先用species进行分组，然后进行分组统计

iris %>% group_by(Species) %>% summarise(avg=mean(Sepal.Length)) %>% arrange(avg) #先用species进行分组，然后进行分组统计,最后依据统计结果进行排序，这样一行代码完成大好多工作，比较搞笑

dplyr::mutate(iris,new=Sepal.Length+Petal.Length)  #mutate函数可以添加新的一列，新列new，是花萼与花瓣长度的总和

a <- data.frame(x1=c("A","B","C"),x2=c(1,2,3))

b <- data.frame(x1=c("A","B","D"),x2=c(T,F,T))

a

b

dplyr::left_join(a,b,by="x1") #对两个数据框进行整合，左链接left_join

dplyr::right_join(a,b,by="x1") #对两个数据框进行整合，左链接left_join,右链接right_join,左链接是以左边的表为基础，因为b没有C行，所以c显示NA

#右链接显示的结果则以右链接为准，a没有D行，所以显示的是NA

dplyr::full_join(a,b,by="x1")

dplyr::inner_join(a,b,by="x1") #内链接是取X1的交集，全链接是取X1的并集

dplyr::semi_join(a,b,by="x1") #半链接是根据右侧表内容对左侧表进行过滤，其实就是取交集

dplyr::anti_join(a,b,by="x1") #反链接也是对右侧表也就是b表进行操作，是将a与b的补集部分输出出来

##几个集合的合并，集合的运算

##以mtcars数据集作为演示

first <- mtcars[c(1:20),]

second <- mtcars[c(10:30),]

first

second

first <- slice(mtcars,1:20)

first <- mutate(first,slice(mtcars,1:20))

first <- mutate(slice(mtcars,1),first)

slice(mtcars,1)

first

first <- slice(mtcars,1:20)

first

mtcars

first1 <- slice(mtcars,1:20)

first1

#以前slice函数可能不会将行名slice下来，可以用mutate函数给first1添加一行

first2 <- mutate(first1,model=rownames(mtcars))

first2 <- mutate(first1,model=rownames(mtcars[1:20],))

first2 <- mutate(first1,model=rownames(mtcars[c(1:20)],))

mtcars[1:20,]

rownames(mtcars[1:20,])

first2 <- mutate(first,model=rownames(mtcars[1:20,]))

first2

first2 <- mutate(model=rownames(mtcars[1:20,]),first1)

first2

first <- slice(mtcars,1:20)

first

second <- slice(mtcars,10:30)

second

intersect(first,second) #intersect函数取交集

dplyr::union_all(first,second) #union_all函数取并集

dplyr::union(first,second) #union函数取非冗余的并集

setdiff(first,second) #取first相对于second的补集

setdiff(second,first) #取second相对于first的补集

##R函数

#R中诸多函数，要学会R至少要掌握400个函数，其中lm函数很好用，可以进行强大的回归分析

#演示，对state.x77进行回归分析，看各州犯罪率与其他条件有无关系

state.x77

stx <- as.data.frame(state.x77) #因为lm操作对象是数据框，所以要先将其转换成数据框再进行操作

stx

calss(stx

)

class(stx)

fit <- lm(Murder~population+Illiteracy+Income+Forst,data=stx)

fit <- lm(Murder~population+Illiteracy+Income+Frost,data=stx)

fit <- lm(Murder~Population+Illiteracy+Income+Frost,data=stx)

fit

summary(fit)

head(stx,10)

##介绍R的返回值

ls() #ls

ls()  #ls返回当前的对象

Sys.Date() #返回当前的系统时间

rm(df)

rm(age) #rm直接删除指定变量，如果成功，没有任何返回

#R中大量的绘图函数直接输出图形，需要单独指定绘图设备，将结果指定到绘图设备中

df

a

plot(a)

age

ex

sdata

plot(sdata)

#plot一下即可将图像绘在右侧白板上，绘图函数对输入数据有要求，比如heatmap函数输入的必须是矩阵，如果是数据框就不行

heatmap(iris)

#lm则必须要数据框，plot函数输入数据类型不同，图形也不同

##选项参数

ls("package:base") #ls函数查看基础base包，可见内有1200多个函数，注意用一个冒号即可

#一般函数的选项参数分为三个部分：输入控制部分，输出控制部分，调节部分；输入部分：file-接一个文件；data-一般指要输入一个数据框；x-表示单独的一个对象，一般都是变量，也可以是矩阵或者列表；x和y-函数需要输入两个变量；x,y,z-函数需要输入3个变量；formula-公式；na.rm-删除缺失值

#比如scatterplot3d是绘制三维图的，所以需要xyz三个变量；

#选项接受哪些参数：main：字符串，不能是向量；na.rm TRUE或者FALSE；axis-side参数只能是1到4；fig-包含4个元素的向量

##数学统计函数--概率函数

rnorm((n=100,mean=15,sd=2) #生成符合要求的一百个随机数

rnorm(n=100,mean=15,sd=2)

round(rnorm(n=100,mean=15,sd=2) ) #使用round函数可以用来取整

x <- rnorm(n=100,mean=15,sd=2)

x

qqnorm(x)

qqnorm(x) #使用qqnorm绘制正态分布图

#R中如何生成随机数

runif #最简单生成随机数的函数，会生成一个0-1之间的随机数

runif(1)

runif(2)

runif(3) #随机生成3个0-1之间的数

q()

#R中有很多数学函数

x <- 1:100

sum(x)

mean(x)

var(x) #var函数是求方差的

sd(x) #sd

sd(x) #sd函数是求标准差的

log(x)

abs(x) #abs函数是求绝对值的

#使用apply系列函数就可以将这些函数应用于矩阵、或者数据框的行和列等进行计算

Sys.Date()

savehistory("~/Desktop/R-data/2021-06-17.Rhistory")

setdiff(first,second) #取first相对于second的补集
setdiff(second,first) #取second相对于first的补集
##R函数
#R中诸多函数，要学会R至少要掌握400个函数，其中lm函数很好用，可以进行强大的回归分析
#演示，对state.x77进行回归分析，看各州犯罪率与其他条件有无关系
state.x77
stx <- as.data.frame(state.x77) #因为lm操作对象是数据框，所以要先将其转换成数据框再进行操作
stx
calss(stx
)
class(stx)
fit <- lm(Murder~Population+Illiteracy+Income+Frost,data=stx)
fit
summary(fit)
head(stx,10)
##介绍R的返回值
ls()  #ls返回当前的对象
Sys.Date() #返回当前的系统时间
rm(df)
rm(age) #rm直接删除指定变量，如果成功，没有任何返回
#R中大量的绘图函数直接输出图形，需要单独指定绘图设备，将结果指定到绘图设备中
df
a
plot(a)
age
ex
sdata
plot(sdata)
#plot一下即可将图像绘在右侧白板上，绘图函数对输入数据有要求，比如heatmap函数输入的必须是矩阵，如果是数据框就不行
heatmap(iris)
#lm则必须要数据框，plot函数输入数据类型不同，图形也不同
##选项参数
ls(package="base") #用ls函数查看base包中的函数
ls(”package=base") #用ls函数查看base包中的函数
ls("package=base") #用ls函数查看base包中的函数
ls("package:base")
ls("package::base") #ls函数查看基础base包，可见内有1200多个函数
ls("package:base") #ls函数查看基础base包，可见内有1200多个函数，注意用一个冒号即可
#一般函数的选项参数分为三个部分：输入控制部分，输出控制部分，调节部分；输入部分：file-接一个文件；data-一般指要输入一个数据框；x-表示单独的一个对象，一般都是变量，也可以是矩阵或者列表；x和y-函数需要输入两个变量；x,y,z-函数需要输入3个变量；formula-公式；na.rm-删除缺失值
#比如scatterplot3d是绘制三维图的，所以需要xyz三个变量；
#选项接受哪些参数：main：字符串，不能是向量；na.rm TRUE或者FALSE；axis-side参数只能是1到4；fig-包含4个元素的向量
##数学统计函数--概率函数
rnorm((n=100,mean=15,sd=2) #生成符合要求的一百个随机数
rnorm(n=100,mean=15,sd=2)
round(rnorm(n=100,mean=15,sd=2) ) #使用round函数可以用来取整
x <- rnorm(n=100,mean=15,sd=2)
x
qqnorm(x)
qqnorm(x) #使用qqnorm绘制正态分布图
#R中如何生成随机数
runif #最简单生成随机数的函数，会生成一个0-1之间的随机数
runif(1)
runif(2)
runif(3) #随机生成3个0-1之间的数
q()
#R中有很多数学函数
x <- 1:100
sum(x)
mean(x)
var(x) #var函数是求方差的
sd(x) #sd函数是求标准差的
log(x)
abs(x) #abs函数是求绝对值的
#使用apply系列函数就可以将这些函数应用于矩阵、或者数据框的行和列等进行计算
Sys.Date()
savehistory("~/Desktop/R-data/2021-06-17.Rhistory")
#选项参数中的三个点表示参数可传递，表示这部分函数的选项参数与其他函数的选项参数是可重合的；有些函数三个点表示没有数量限制，比如data.frame表示可以将多个向量进行合并，没有数量限制；另有函数包含formula参数，表示输入内容应该写成一个公式，比如加减乘除，一般是波浪线连接的公式，波浪线表示相关；输入还有na.rm和na.omit，如果想处理缺失值还需要这两个参数，这些是主要的输入选项参数
#main参数是用来设置标题的，只能赋字符串
#R中常用的数学统计函数
rnorm(n=100,mean=15,sd=2) #rnorm随机数函数，生成均值15，标准差为2的100个随机数
round(rnorm(n=100,mean=15,sd=2)) #round函数对随机数进行取整
x <- rnorm(n=100,mean=15,sd=2) #得到正态分布数据
qqnorm(x) #使用qqnorm绘制正态分布图
x <- rnorm(n=100,mean=15,sd=2) #得到正态分布数据
x
qqnorm(x)
runif(1) #生成一个0-1之间的数
runif(2) #生成2个0-1之间的数
runif(10)*10 #随机生成10个0-10之间的数，也就是在原来随机数的基础上*10，但是这样并不高效，还有其他解法
runif(50,min = 1,max = 100)  #runif函数内置min和max参数；该代码表示生成50个1-100之间的随机数
dgama(c(1:9),sahpe)
dgama(c(1:9),sahpe=2,rate=1) #随机生成的概率密度数字
?sgama
?dgama
set.seed(666)
runif(4)
set.seed(666)
runif(4)
##set.seed函数括号里面加一个数字，可以随便添加一个数字，runif得到随机数之后，再运行一次set.seed函数，再运行runif，就可以得到和之前相同的随机数，避免了每次得到的随机数不同的问题
###描述性统计函数
myvars <- mtcars[c("mpg","hp","wt","am")]
summary(myvars) #summary函数运行一次就可以对这个数据集进行详细统计
?fivenum
fivenum(myvars$hp) #fivenum返回的5个数是最小值、四分位低值，中位数、四分位高值、最大值
library(Hmisc)
install.packages("Hmisc")
library(Hmisc)
?describe
describe(myvars) #返回这四列数据的前五个最大最小值，中位数，以及四分位数和四分之三位数，缺失值的个数，行数（Observation)，列数(Variables)
myvars
install.packages("pastecs")
library(pastecs)
stat.desc(myvars) #pastecs包中的stat.desc函数可以进行多种描述性统计函数
stat.desc(myvars,basic=TRUE) #some basic result
##nbr.val means the number of all result; nrb.null means the number of null ;nbr.na means the nubmer of NA;
stat.desc(myvars,nrom=TRUE) #some norm result
stat.desc(myvars,norm = TRUE) #some norm result
library(psych)
install.packages("psych")
library(psych)
describe(myvars) #describe function can provide number of no-NA
describe(myvars) #describe function can provide number of no-NA and other results
describe(myvars) #describe function can provide number of no-NA and other results such as ending mean result ,this result means number of mean result without the max and min numerical value
describe(myvars,trim = 0.1) # trim parameter means to remove the max and min 10% part
# there are describe function in both Hmisc and psych  package , the way to determine which function is using now is that R consider the latter one is using now
Hmisc::describe(myvars) # the way to use the previous function other than the latter one
library(MASS)
Cars93
#Cars93 is a data set including index of many cars manufactured in 1993
head(Cars93)
aggregate(Cars93[c("Min.Price","Price","Max.Price","MPG.city")],by = list(Manufacturer=Cars93$Manufacturer),mean)  #calculate mean number of Min.Price etc using Manufacturer as grouping standard , using "by" parameter
#note the parameters of  aggregate function
aggregate(Cars93[c("Min.Price","Price","Max.Price","MPG.city")],by = list(Manufacturer=Cars93$Manufacturer),sum)
aggregate(Cars93[c("Min.Price","Price","Max.Price","MPG.city")],by = list(Orign=Cars93$Orign,Manufacturer=Cars93$Manufacturer),mean)  #calculate mean number of Min.Price etc using Manufacturer as grouping standard , using "by" parameter
aggregate(Cars93[c("Min.Price","Price","Max.Price","MPG.city")],by = list(Origin=Cars93$Origin,Manufacturer=Cars93$Manufacturer),mean)  #calculate mean number of Min.Price etc using Manufacturer as grouping standard , using "by" parameter
aggregate(Cars93[c("Min.Price","Price","Max.Price","MPG.city")],by = list(Origin=Cars93$Origin,Manufacturer=Cars93$Manufacturer),mean)  #calculate mean number of Min.Price etc using Manufacturer and Origin as grouping standards , using "by" parameter
#the defect of aggregate function is that only one function such as sum can be used once
library(doBy)
install.packages("doBy")
library(doBy)
summaryBy(mpg+hp+wt ~ am,data=myvars, FUN = mean) #mpg+hp+wt : use only col names; ~ means correlation ;am means grouping stardards ; FUN means function ;
summaryBy(mpg+hp+wt ~ am,data=myvars, FUN = mean+sum) #mpg+hp+wt : use only col names; ~ means correlation ;am means grouping stardards ; FUN means function ;
summaryBy(mpg+hp+wt ~ am,data=myvars, FUN = c(mean+sum) #mpg+hp+wt : use only col names; ~ means correlation ;am means grouping stardards ; FUN means function ;
)
summaryBy(mpg+hp+wt ~ am,data=myvars, FUN = c(mean,sum) #mpg+hp+wt : use only col names; ~ means correlation ;am means grouping stardards ; FUN means function ;
)
# use two function such as mean and sum simultaneouly usinf c();
psych::describe.by(myvars,list(am=myvars$am)) # use describe function belonged to psych package ,list means gouping standards, can get many kinds of statistics
psych::describe.by(myvars,list(am=myvars$am,manufactuer=myvars$Manufacturer))
psych::describe.by(myvars,list(am=myvars$am,Manufacturer=myvars$Manufacturer))
##Frequency statistics Function
mtcars
cyl <- as.factor(mtcars$cyl) #using cyl col as factor to count frequency , transfer cyl col to factor
cyl
split(mtcars,cyl) # using split function to group the mtcars (cyl as group standards)
cut(mtcars$mpg,c(10,50,10)) #using cut function to cut the mpg col ,cut 10-50 to four parts , standards=10;
cut(mtcars$mpg,c(seq(10,50,10)) #using cut function to cut the mpg col ,cut 10-50 to four parts , standards=10;
)
table(mtcars$cyl) #using table function to do frequency statistics
table(cut(mtcars$mpg,c(seq(10,50,10))))
table(cut(mtcars$mpg,c(seq(10,50,10)))) # using table function to do frequency statistics of mtcars$mpg , cut mpg into four parts previous
table(cut(mtcars$mpg,c(seq(10,50,10)))) # using table function to do frequency statistics of mtcars$mpg , cut mpg into four parts previously
prop.table(table(mtcars$cyl)) # using prop.table function to do proportion statistics
prop.table(table(cut(mtcars$mpg,c(seq(10,50,10)))))
#how to deal with 2 dimensions data.frame and get two dimension frequency
library(vcd)  #using Arthritis data set to example the function
install.packages("grid")
library(vcd)
Arthritis # load the Arthritis dataset
head(Arthritis)
table(Arthritis$Treatment,Arthritis$Improved) # count the frequency of this dataset
# treatment including Placebo and Treated , improved including None,some,Marked as two factor
# get a contingency table
with(data = Arthritis,(table(Treatment,Improved)))  # key to deal with too many variables  # do not forget the , after data = Arthritis
?xtabs
xtabs(~Treatment+Improved,data = Arthritis)  # use formula in this xtab function , the first is formula ,~T+I means use these two factor as statistical indicators
?margin.table
x <- xtabs(~Treatment+Improved,data = Arthritis)
margin.table(x)
margin.table(x,1)
margin.table(x,2)
margin.table(x,2) # margin.table() , first parameter means the object,second parameter : 1 means row ; 2 means col;
#margin.table() function means analyze the two dimensions contingency table according to only rows or cols ; not two of them
x
prob.table(x,1)
prop.table(x,1)
prop.table(x,2)  # use prop.table function to get percentage of different events accrording to rows or cols
x
x
addmargins(x)
addmargins(x) # use addmargins function to get sum of rows or cols,1 means row,2 means col
addmargins(x,1)
addmargins(x,2)
xtabs(~Treatment+Improved+Sex,data = Arthritis)  # creat a three dimensional contingency table
y <- xtabs(~Treatment+Improved+Sex,data = Arthritis)
ftable(y)  # use ftable function to convert the three dimensional contingency table to a flat contingency table
##Independency test
# three independency test algorithm : Chi-square test; Fisher test ;Cochran-Mantel-Haenszel test
# example: get if it is independent between two variables: Treatment , Improved
mytable <- xtabs(~Treatment+Improved,data = Arthritis) # get the frequency
mytable
?chisq.test
chisq.test(mytable)  # use chisq.test() function to do Chi-Square test
mytable <- xtabs(~Sex+Improved,data = Arthritis)
chisq.test(mytable)
#mytable <- xtabs(~Sex+Improved,data = Arthritis)  # no order difference between the two varibales:Sex/Improved
?fisher.test
## in fisher exact test , row and col are independent which is different from chisq.test
fisher.test(mytable)
mytable
#mytable <- xtabs(~Treatment+Improved,data = Arthritis)  # no order difference between the two varibales:Sex/Improved
fisher.test(mytable)
mytable <- xtabs(~Treatment+Improved,data = Arthritis)
mytable
fisher.test(mytable)
mantelhaen.test(mytable) #两个名义变量在第三个变量每一层中都是条件独立的
mytable <- xtabs(~Treatment+Improved+Sex,data = Arthritis)
mantelhaen.test(mytable)
## order has influence in mantelhaen.test() function
##this mytable means sex is the grouping standard ; result means treatment and improved are not independent in every levels of sex (male and female )
##correlation analysis functions
## correlation coefficient : Pearson correlation coefficient;Spearman correlation coefficient ; Kendall correlation coefficient;polychoric correlation coefficient;polyserial correlation coefficient; Partial correlation coefficient;
##use cor function to count three kinds of correlation coefficient
?cor
##default is the Pearson correlation coefficient
## default is the Pearson correlation coefficient
state.x77
head(state.x77)
class(state.x77)
cor(state.x77)  # get correlation coefficient between two variables; such as the correlation coefficient between income and crime in Amlabama
## too easy to get the correlation coefficient
cor(state.x77)
head(state.x77)
## cov() function to get the covariation between 2 variables ; means the total error of two variables
## when count Partial correlation coefficient , Covariation result is needed
cov(state.x77)
## cov() and cor() functions influence the same meaning
## next:how to analyze relationship between one variable and the other variable
colnames(state.x77)
x <- state.x77[,c(1,2,3,6)]
x
head(x)
y <- state.x77[,c(4,5)]
head(y)
cor(x,y)  # result including correlationship is easy and clean
library(ggm)
install.packages("ggm")
library(ggm)
## install and load the expansion pack of R :ggm , which can be used to count Partial correlation coefficient
## partial correlation means when one or more variables are controlled , the correlation coefficient between the remaining variables
pcor(c(1,5,2,3,6),cov(state.x77))
?pcor()
cov(state.x77)
## how to do significant analytical test when getting the partial correlation coefficient
、cor.test
?cor.test
cor.test(state.x77[,3],state.x77[,5])
colnames(state.x77)
cor.test(state.x77[,2],state.x77[,5])
## confidence interval
## cor.test can only analyze one group of variables , corr.test in psych package can analyze two more variables once
?corr.test
corr.test(state.x77)
corr.test(state.x77,short = FALSE)
corr.test(state.x77)
pcor.test(xtate.x77)
pcor.test(state.x77)
x <- pcor(c(1,5,2,3,6),cov(state.x77))
x
pcor.test(x,3,,50)
pcor.test(x,3,50)
## correlation test of packet data
libarary(MASS)
library(MASS)
UScrime
head(UScrime)
?t.test
t.test(Prob~So,data = UScrime)
t.test(Prob~So,data = UScrime) # example of two teams
###2021-06018 Drawing functions
## cor.test can only analyze one group of variables , corr.test in psych package can analyze two more variables once
## many drawing functions in R and its expansion package : basic drawing system ; lattice package;ggplot2 package;grid package
ls("package::graphcis")
ls("package::graphics")
ls("package:graphics")
## get drawing function in graphics pack which is default loaded when R stary
head(state.x77)
class(state.x77)
heatmap(state.x77)
hist(state.x77)
stars(state.x77)
demo(graphics) # show some plot that R can draw
## the most common functions : plot and par
head(women)
?plot
plot(women$height)
plot(women$height,women$weight)
plot(as.factor(women$height))  # get histogram when object is factor
plot(as.factor(women$weight))  # get histogram when object is factor
plot(mtcars$cyl)
plot(as.factor(mtcars$cyl))
head(mtcars)
plot(as.factor(mtcars$qsec))
plot(as.factor(mtcars$disp))
plot(as.factor(mtcars$cyl),mtcars$carb)
class(mtcars$carb)
plot(as.factor(mtcars$cyl),mtcars$gear)
plot(mtcars$carb,as.factor(mtcars$cyl))
plot(as.factor(mtcars$carb),as.factor(mtcars$cyl))
plot(as.factor(mtcars$cyl,as.factor(mtcars$carb)))
plot(as.factor(mtcars$cyl),as.factor(mtcars$carb))
plot(women$height~weight) # get the relationship of two variables: height and weight
plot(women$height~women$weight) # get the relationship of two variables: height and weight
fit <- lm(women$height,women$weight)
fit <- lm(height~weight,data = women)
fit
plot(fit)
plot(fit)
method(plot)
Method(plot)
methods(plot)
?par
### par function is used to set up the drawing parameter
par()
plot(as.factor(mtcars$cyl),col=c("red","green","blue"))
plot(as.factor(mtcars$cyl),col=c("pink","green","blue"))
plot(as.factor(mtcars$cyl),col=c("pink","green","yellow"))
plot(as.factor(mtcars$cyl),col=c("pink","gray","yellow"))
plot(as.factor(mtcars$cyl),col=c("pink","gray","light yellow"))
## custom function
cor
mtcas
head(mtcars)
mtcars$mpg
split(mtcars$mpg,10,30,5)
?split
?data.table:split
mystats <- function(x,na.omit=FALSE){
if(na.omit)
x <- x[!is.na(x)]
m <- mean(x)
n <- length(x)
s <- sd(x)
skew <- sum((x-m^3)-s^3)/n
kurt <- sum((x-m^4)-s^4)/n-3
return(c(n=n,mean=m,std=s,skew=skew,kurtosis=kurt))
}
mystats
x <- 1:100
mystats(x)
for (i in 1:10){print("hello,world")}
i = 1
i=1,while(i <=10),{print ("hello,world")}
q()
women
plot(women)
fit <- lm(weight~height,data = women)
fit
summary(fit)
summary.lm(fit)
q()
fit(women)
librbary(ggm)
library(ggm)
women
fit(women)
library(MASS)
fit
fit(women)
plot(women)
lm(women$height,women$weight)
ls("package:graphics")
fit=lm(women$height,women$weight)
women
plot(women)
plot(women$height,women$weight)
?lm
fit <- lm(weight~height,data = women)
fit
coefficients(fit)
confint(fit) ## get the confidence interval
confint(fit,level = 0.5) ## get the confidence interval, adjust it using the "level" options
confint(fit,level = 0.25) ## get the confidence interval, adjust it using the "level" options
fitted(fit)  #get the predictive value of the fitting model
women$weight-fitted(fit) # get the residual using this formula
residuals(fit) # can also get the residual value
women1 <- women
predict(fit,women1) # predcit the new dataset using our fitting model
women1
plot(fit) #get the image of this fitting model using plot function
# we get four images to describe this fitting model
plot(women$height,women$weight)
abline(fit)
fit2 <- lm(weight~height+I(height^2),data = women)  # if we want to get height2 in this function ,it should be I(height^2),bacause ^ in this formula have other meanings
fit2
summary(fit2)
plot(women$height,women$weight)
abline(fit) # get the result of the first linear regression
lines(women$height,fitted(fit2),col = "light red") # get the result of the second regression which is might not be straight line
lines(women$height,fitted(fit2),col = "red") # get the result of the second regression which is might not be straight line
lines(women$height,fitted(fit2),col = "red") # get the binomial regression model which is more accurate
fit3 <- lm(weight~height+I(height^3),data = women)
lines(women$height,fitted(fit3),col = "yello")
lines(women$height,fitted(fit3),col = "yellow")
fit3 <- lm(weight~height+I(height^2)+I(height^3),data = women)
lines(women$height,fitted(fit3),col = "yellow")
summary(fit3)
###2021-06-19 multiple linear regression
head(state.x77)
states <- as.data.frame(state.x77)
class(states)
# transfer state.x77 to data.frame because it must be format of data.frame when using lm function
states <- as.data.frame(state.x77[,c(5,1,3,2,7)])
states
head(states)
fit <- lm(Murder~Population+ Illiteracy+ Income+ Frost,data = states)
summary(fit)
coef(fit)
head(mtcars)
fit <- lm(mpg~hp+ wt+hp:wt,data = mtcars) #hp:wt means there is a relationship between hp and wt which we still dont know
summary(fit)
## if we have four variables , there will be many predictive models , here shows how to choose the most accurate model
?AIC
fit1 <- lm(Murder~Population+ Illiteracy+ Income+ Frost,data = states)
fit2 <- lm(Murder~Population+ Illiteracy,data = states)
AIC(fit1,fit2)
## the AIC result shows fit2 model is more accurate cause the AIC value 241.6429\237.6565 , fit2' value 237.6565 is smaller;
## if there are more models/ variables , it is no efficiency to use AIC function , now we can use stepwise regression and full-subset regression
library(MASS)
?stepAIC #use stepAIC function in MASS package to do stepwise `regular expression`
install.packages("leaps")
library(leaps)
?regsubsets # use regsubsets function in leaps pack to do full-subsets regression
## Regression diagnosis
## plot(result of lm function) to create 4 images to evaluate fitting model
fit <- lm(weight~height,data=women)
plot(fit)  # create images to evaluate fitting model
par(mfrow=c(2,2))  # set up the dsitribution of these 4 images,2*2
plot(fit)
## how to eatimate the independence of the fitting model
fit2 <- lm(weight~height+I(height^2),data=women)
opar <- par(no.readonly = TRUE)
par(mfrow=c(2,2))
plot(fit2)
### Analysis of Variance 简称ANOVA
library(multcomp)  # 50 cholesterol patients treated by different medicines
install.packages("multcomp")
library(multcomp)
cholesterol
attach(cholesterol)
table(trt)
?aggregate
aggregate(response,by=list(trt),FUN=mean())
aggregate(response,by=list(trt),FUN=mean)  # now we get the mean value of the response value grouping by treating times:1 time;2 times;4times and other optional drug D and E
fit <- aov(response~ trt,data = cholesterol) #using aov function to do the variance analysis of this dataset/experimental result
fit
summary(fit)
## F value bigger ; difference between teams is more significant;
plot(fit)
fit.lm <- aov(response~trt,data = choleaterol)  # use lm
fit.lm <- lm(response~trt,data = choleaterol)  # use lm function to do variance analysis , and analyze the difference between aov and lm functions
fit.lm <- lm(response~trt,data = cholesterol)  # use lm function to do variance analysis , and analyze the difference between aov and lm functions
fit.lm
summary(fit.lm)
head(litter)
table(litter$dose)
attach(litter)  # the right code using aggregate function
aggreagate(weight,by=list(dose),FUN=mean)
aggregate(weight,by=list(dose),FUN=mean)
# x mean value between teams is different ; use variance analysis to test if there is significant difference different team
?aov
fit <- aov(weight~gesttime+dose,data = litter) # use aov function to do variance analysis;gesttime is vovariate and dose is independent variables
summary(fit)
# this result indicates that gesttime and dose influences the birth weight ; gesttime is more important and dose also affects ; see the F value and P value (p value is the eatimate stardand to estimate F value ); also see the **
## Two way analysis of variance - Example
savehistory("~/Desktop/R-data/2021-06-18.Rhistory")

2021-06-24

> # scatter plot
>  getwd()
[1] "/Users/yangyi"
> m <- read.table("Rdata/Rplotdata/prok_representative.txt"，sep = "\t")
Error: unexpected input in "m <- read.table("Rdata/Rplotdata/prok_representative.txt"，"
> m <- read.table("Rdata/Rplotdata/prok_representative.txt",sep = "\t")
Error in file(file, "rt") : cannot open the connection
In addition: Warning message:
In file(file, "rt") :
  cannot open file 'Rdata/Rplotdata/prok_representative.txt': No such file or directory
> getwd()
[1] "/Users/yangyi"
> setwd("/Users/yangyi/Desktop/Rdata")
> getwd()
[1] "/Users/yangyi/Desktop/Rdata"
> m <- read.table("Rplotdata/prok_representative.txt",sep = "\t")
> view()
Error in view() : could not find function "view"
> view(m)
Error in view(m) : could not find function "view"
> View(m)
> genome_size <- m[,2]
> gene_size <- m[,4]
> head(genome_size)
[1] 2.848380 7.477260 4.557580 5.967560 0.159662 5.369770
> head(gene_size)
[1] 2823 6618 3904 5311  213 4812
> plot(genome_size,gene_size )
> plot(genome_size,gene_size,pch=16 )
> plot(genome_size,gene_size,pch=16 ) # pch parameter sets the shape of the dot
> plot(genome_size,gene_size,pch=16,cex )
Error in plot.xy(xy, type, ...) : object 'cex' not found
> plot(genome_size,gene_size,pch=16,cex=0.8 )  # cex parameter sets the size of the dot
> fit <- lm(genome_size,gene_size)
Error in formula.default(object, env = baseenv()) : invalid formula
> fit <- lm(genome_size~gene_size)
> fit

Call:
lm(formula = genome_size ~ gene_size)

Coefficients:
(Intercept)    gene_size  
  -0.203979     0.001147  

> fit <- lm(gene_size~genome_size)
> fit

Call:
lm(formula = gene_size ~ genome_size)

Coefficients:
(Intercept)  genome_size  
      286.6        843.7  

> abline(fit,col="blue",lwd=1.8)  # add the regression curve using abline function
> summary(fit)

Call:
lm(formula = gene_size ~ genome_size)

Residuals:
     Min       1Q   Median       3Q      Max
-10606.7   -102.1     18.0    116.8   2551.7

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  286.598     11.211   25.57   <2e-16 ***
genome_size  843.691      2.608  323.55   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 294.4 on 3502 degrees of freedom
Multiple R-squared:  0.9676,    Adjusted R-squared:  0.9676
F-statistic: 1.047e+05 on 1 and 3502 DF,  p-value: < 2.2e-16

> ## the formula is y=843.691x+286.598   R^=0.9676
> text(12,10000,labels="y=843.691+286.598\nR2=0.97")  # add the fromula in the picture , 12 means x coordinate,10000 means y coordinate, nR2 means R^2
> text(12,10000,labels="y=843.691+286.598\nR2=0.97",col="red")
> summary(fit)$adj.r.square  # get the R^2 value from fit result
[1] 0.9676204
> rr <- round(summary(fit)$adj.r.square,2)  # keep two significant digits
> rr
[1] 0.97
> intercept <- summary(fit)$Coefficients[1,2]
> intercept
NULL
> intercept <- round(summary(fit)$Coefficients[1],2) # get the intercept of the fit summary
Error in round(summary(fit)$Coefficients[1], 2) :
  non-numeric argument to mathematical function
> intercept <- round(summary(fit)$coefficients[1],2) # get the intercept of the fit summary
> intercept <- round(summary(fit)$coefficients[1],2) # get the intercept of the fit summary using fit$coefficients function
> intercept
[1] 286.6
> slope <- round(summary(fit)$coefficients[2],2) # get theslope of the fit summary using fit$coefficients[2] function
> slope
[1] 843.69
> ## then compose the formula using rr,intercept,slope
> eq <- bquote(atop("y="*(slope)*"x+"*(intercept),R^2==.(rr)))  # use atop and bquote function to create the function
> eq
atop("y=" * (slope) * "x+" * (intercept), R^2 == 0.97)
> text(12,6e3,eq)
> eq <- bquote(atop("y="*.(slope)*"x+"*.(intercept),R^2==.(rr)))  # use atop and bquote function to create the function
> text(12,6e3,eq)
> plot(genome_size,gene_size,pch=16,cex )
Error in plot.xy(xy, type, ...) : object 'cex' not found
> plot(genome_size,gene_size,pch=16,cex=0.8 )
> fit <- lm(gene_size~genome~)
Error: unexpected ')' in "fit <- lm(gene_size~genome~)"
> fit <- lm(gene_size~genome_size)
> summary(fit)