R语言学习笔记-三

R for Data Science

R语言实战 学习笔记

chapter 1 R语言介绍

1.1 帮助函数

函数 功能
help.start() 打开帮助文档首页
help("foo")或?foo 查看函数foo的帮助（引号可以省略）
help.search("foo")或??foo 以foo未关键词搜索本地帮助文档
example("foo") 函数foo的使用实例（引号可以省略）
RSiteSearch("foo") 以foo为关键词搜索在线文档和邮件列表存档
apropos("foo",mode="function") 列出名称中含有foo的所有可用函数
data() 列出当前已加载包中所有可用示例数据集
vignette() 列出当前已安装包中所有可用的vignette文档
vignette("foo") 为主题foo显示指定的vignette文档

1.2 管理R工作空间的函数

函数 功能
getwd() 显示当前的工作目录
setwd("mydirectory") 修改当前的工作目录为mydirectory
ls() 列出当前工作空间中的对象
rm(objectlist) 移除（删除）一个或多个对象
help(options) 显示可用选项的说明
options() 显示或设置当前选项
history(#) 显示最近使用过的#个命令（默认值为25）
savehistory("myfile") 保存命令历史到文件myfile中（默认值为.Rhistory）
loadhistory("myfile") 载入一个命令历史文件（默认值为.Rhistory）
save.image("myfile") 保存工作空间到文件myfile中（默认值为.RData）
save(objectlist, file="myfile") 保存指定对象到一个文件中
load("myfile") 读取一个工作空间到当前会话中（默认值为.RData）
q() 退出R。将会询问你是否保存工作空间

1.3 输入和输出

输入
source("filename")可在当前会话中执行一个脚本
文本输出
sink("filename")将输出重定向到文件filename，append=TRUE，可追加；split=TRUE可将输出同时发送到屏幕和输出文件中；不加参数调用仅向屏幕返回输出结果。
图形输出

函数 功能
bmp("filename.bmp") BMP文件
jpeg("filename.jpg") JPEG文件
pdf("filename.pdf") PDF文件
png("filename.png") PNG文件
postscript("filename.ps") PostScript文件
svg("filename.svg") SVG文件
win.metafile("filename.wmf") Windows图元文件

1.4 包

R语言包下载

# 安装包
install.packages()
# 更新包
update.packages()
# 列出已安装包，版本号及依赖关系等
installed.packages()
# 载入包
library()

# 通过车身重量（wt）预测每加仑行驶的英里数（mpg）线性回归
lm(mpg~wt, data=mtcars)
# 保存为对象
lmfit <- lm(mpg~wt, data=mtcars)
# 显示分析结果的统计概要
summary(lmfit)
# 生成回归诊断图形
plot(lmfit)

chapter 2 创建数据集

# 创建一个5*4的矩阵
y <- matrix(1:20, nrow=5, ncol=4)

cells <- c(1,26,24,68)
rnames <- c("R1", "R2")
cnames <- c("C1", "C2")
# 按行填充2*2矩阵
mymatrix <- matrix(cells, nrow=2, ncol=2, byrow=TRUE, dimnames=list(rnames, cnames))

# 按列填充2*2矩阵
mymatrix <- matrix(cells, nrow=2, ncol=2, byrow=FALSE, dimnames=list(rnames, cnames))

# 创建一个数组  myarray <- array(vector, dimensions, dimnames)
dim1 <- c("A1", "A2") 
dim2 <- c("B1", "B2", "B3") 
dim3 <- c("C1", "C2", "C3", "C4")
z <- array(1:24, c(2, 3, 4), dimnames=list(dim1, dim2, dim3))

# 创建一个数据框 mydata <- data.frame(col1, col2, col3,...)
patientID <- c(1, 2, 3, 4)
age <- c(25, 34, 28, 52)
diabetes <- c("Type1", "Type2", "Type1", "Type1")
status <- c("Poor", "Improved", "Excellent", "Poor")
patientdata <- data.frame(patientID, age, diabetes, status)

# 取列数据
patientdata[1:2]
patientdata[c("diabetes","status")]

# 取一列数据
patientdata$age

# 以patientdata$diabetes为行，以patientdata$status为列，作一个列联表
table(patientdata$diabetes, patientdata$status)

# 获取每加仑行驶英里数（mpg）变量的描述性统计量
summary(mtcars$mpg)
# 绘制此变量与发动机排量（disp）的散点图
plot(mtcars$mpg, mtcars$disp)
# 绘制此变量与车身重量（wt）的散点图
plot(mtcars$mpg, mtcars$wt)

# attach()可将数据框添加到R的搜索路径中
# detach()将数据框从搜索路径中移除
# 名称相同时存在冲突，原始对象取得优先权
# 上述三句代码可简化为
attach(mtcars)
 summary(mpg)
 plot(mpg, disp)
 plot(mpg, wt)
detach(mtcars)

# 使用with()语法进行平替，使用<-产生局部变量；<<-产生全局变量
with(mtcars, {
 print(summary(mpg))
 plot(mpg, disp)
 plot(mpg, wt)
})

变量可归结为名义型、有序型或连续型变量。

  • 名义型变量是没有顺序之分的类别变量，eg，patientdata$diabetes。
  • 有序型变量表示一种顺序关系，而非数量关系，eg，patientdata$status。
  • 连续型变量可以呈现为某个范围内的任意值，并同时表示了顺序和数量，eg，patientdata$age。

类别（名义型）变量和有序类别（有序型）变量在R中称为因子（factor）

diabetes <- c("Type1", "Type2", "Type1", "Type1")
status <- c("Poor", "Improved", "Excellent", "Poor")
# 函数factor()以一个整数向量的形式存储类别值，整数的取值范围是[1...k]
# ordered=TRUE <=> 因子的水平默认依字母顺序创建
status <- factor(status, ordered=TRUE)
# 指定levels选项来覆盖默认排序
status <- factor(status, order=TRUE, levels=c("Poor", "Improved", "Excellent"))

# 创建一个列表
g <- "My First List"
h <- c(25, 26, 18, 39)
j <- matrix(1:10, nrow=5)
k <- c("one", "two", "three")
mylist <- list(title=g, ages=h, j, k)

# mylist[2] 和 mylist[[2]] 获取的不一样

# 创建含有三个变量的mydata数据框
mydata <- data.frame(age=numeric(0), gender=character(0), weight=numeric(0))
# 使用编辑器进行数据输入，或使用fix(mydata)
mydata <- edit(mydata)
# or use code
mydatatxt <- "
age gender weight
25 m 166
30 f 115
18 f 120
"
mydata <- read.table(header=TRUE, text=mydatatxt)

# 从带分隔符的文本文件导入数据 mydataframe <- read.table(file, options)

2.1 函数read.table()的选项

选项 描述
header 一个表示文件是否在第一行包含了变量名的逻辑型变量
sep 分开数据值的分隔符。默认是 sep=""，这表示了一个或多个空格、制表符、换行或回车。使用 sep=","来读取用逗号来分隔行内数据的文件，使用 sep="\t"来读取使用制表符来分割行内数据的文件
row.names 一个用于指定一个或多个行标记符的可选参数
col.names 如果数据文件的第一行不包括变量名（header=FASLE），你可以用 col.names 去指定一个包含变量名的字符向量。如果 header=FALSE 以及 col.names 选项被省略了，变量会被分别命名为 V1、V2，以此类推
na.strings 可选的用于表示缺失值的字符向量。比如说，na.strings=c("-9", "?")把-9 和?值在读取数据的时候转换成 NA
colClasses 可选的分配到每一列的类向量。比如说，colClasses=c("numeric", "numeric", "character", "NULL", "numeric")把前两列读取为数值型变量，把第三列读取为字符型向量，跳过第四列，把第五列读取为数值型向量。如果数据有多余五列，colClasses 的值会被循环。当你在读取大型文本文件的时候，加上 colClasses 选项可以可观地提升处理的速度
quote 用于对有特殊字符的字符串划定界限的字符串。默认值是双引号（"）或单引号（'）
skip 读取数据前跳过的行的数目。这个选项在跳过头注释的时候比较有用
stringAsFactors 一个逻辑变量，标记处字符向量是否需要转化成因子。默认值是 TRUE，除非它被 colClases所覆盖。当你在处理大型文本文件的时候，设置成 stringsAsFactors=FALSE 可以提升处理速度
text 一个指定文字进行处理的字符串。如果 text 被设置了，file 应该被留空。 
# studentgrades.csv
StudentID,First,Last,Math,Science,Social Studies 
011,Bob,Smith,90,80,67 
012,Jane,Weary,75,,80 
010,Dan,"Thornton, III",65,75,70 
040,Mary,"O'Leary",90,95,92

# 以下两种方式都可以读入.csv数据，但是使用第一种也没办法得到factor的数据
grades <- read.table("studentgrades.csv", header=TRUE, row.names="StudentID", sep=",")

grades <- read.table("studentgrades.csv", header=TRUE, row.names="StudentID", sep=",",  colClasses=c("character", "character", "character", "numeric", "numeric", "numeric"))

2.2 使用连接到导入数据

例如，函数file(),gzfile(),bzfile(),xzfile(),unz(),url()

2.3 导入数据

  • 1 导入excel数据

需要xlsx、xlsxjars、rJava包
==>安装存在问题(to be continued)

  • 2 导入XML数据

XML的相关文档

  • 3 从网页抓取数据

readLines()
grep()
gsub()
RCurl包和XML包
Programming With R
twitterR获取Twitter数据
Rfacebook获取Facebook数据
Rflickr获取Flicker数据
Web Technologies and Services全面的列表

  • 4 导入SPSS数据

使用foregign包的read.spss()；使用Hmisc包的spss.get()

install.packages("Hmisc")
library(Hmisc)
mydataframe <- spss.get("mydata.sav", use.value.labels=TRUE)
  • 5 导入SAS数据

使用foregign包的read.ssd()；使用Hmisc包的sas.get()；使用sas7dbat包的read.sas7bdat()；SAS中的sas.get()

# 未验证
library(Hmisc) 
datadir <- "C:/mydata" 
sasexe <- "C:/Program Files/SASHome/SASFoundation/9.4/sas.exe" 
mydata <- sas.get(libraryName=datadir, member="clients", sasprog=sasexe)
  • 6 导入Stata数据
library(foreign)
mydataframe <- read.dta("mydata.dta")
  • 7 导入NetCDF数据
library(ncdf) 
nc <- nc_open("mynetCDFfile") 
myarray <- get.var.ncdf(nc, myvar)
  • 8 导入HDF5数据
source("http://bioconductor.org/biocLite.R") 
biocLite("rhdf5")
  • 9 访问数据库管理系统
    • ODBC接口
      安装RODBC包
    library(RODBC) 
myconn <-odbcConnect("mydsn", uid="Rob", pwd="aardvark") 
crimedat <- sqlFetch(myconn, Crime) 
pundat <- sqlQuery(myconn, "select * from Punishment") 
close(myconn)
    • DBI相关包
      RMySQL\ROracle\RPOSTgreSQL\RSQLite
函数 描述
odbcConnect(dsn,uid="",pwd="") 建立一个到 ODBC 数据库的连接
sqlFetch(channel,sqltable) 读取 ODBC数据库中的某个表到一个数据框中
sqlQuery(channel,query) 向 ODBC 数据库提交一个查询并返回结果
sqlSave(channel,mydf,tablename=sqtable,append=FALSE) 将数据框写入或更新（append=TRUE）到 ODBC数据库的某个表中
sqlDrop(channel,sqtable) 删除 ODBC 数据库中的某个表
close(channel) 关闭连接
  • 10 通过Stat/Transfer导入数据

2.4 数据集的标注

  • 1 变量标签
names(patientdata)[2] <- "Age at hospitalization (in years)"
  • 2 值标签
patientdata$gender <- factor(patientdata$gender, levels = c(1,2), labels = c("male", "female"))

2.5 处理数据对象的使用函数

函数 用途
length(object) 显示对象中元素/成分的数量
dim(object) 显示某个对象的维度
str(object) 显示某个对象的结构
class(object) 显示某个对象的类或类型
mode(object) 显示某个对象的模式
names(object) 显示某对象中各成分的名称
c(object, object,...) 将对象合并入一个向量
cbind(object, object, ...) 按列合并对象
rbind(object, object, ...) 按行合并对象
object 输出某个对象
head(object) 列出某个对象的开始六行
tail(object) 列出某个对象的最后六行
ls() 显示当前的对象列表
rm(object, object, ...) 删除一个或更多个对象。语句 rm(list = ls())将删除当前工作环境中的几乎所有对象
newobject <- edit(object) 编辑对象并另存为 newobject
fix(object) 直接编辑对象

chapter 3 图形初阶

# 交互式风格
attach(mtcars)
plot(wt, mpg)
abline(lm(mpg~wt))
title("Regression of MPG on Weight")
detach(mtcars)

# 将图片保存为pdf
pdf("mygraph.pdf")
 attach(mtcars)
 plot(wt, mpg)
 abline(lm(mpg~wt))
 title("Regression of MPG on Weight")
 detach(mtcars)
dev.off()

# 在创建一幅新图形之前打开一个新的图形窗口-伪代码
dev.new()
statements to create graph 1
dev.new()
statements to create a graph 2
etc.
# 函数dev.new()、dev.next()、dev.prev()、dev.set()和dev.off()同时打开多个图形窗口

# 一个简单的图形例子
dose <- c(20, 30, 40, 45, 60)
drugA <- c(16, 20, 27, 40, 60)
drugB <- c(15, 18, 25, 31, 40)
plot(dose, drugA, type="b")

# 使用par()指定图形参数，除再次修改，会话结束前一直有效，语法：par(optionname=value, optionname=name,...)
# 添加参数no.readonly=TRUE可以生成一个可以修改的当前图形参数列表
opar <- par(no.readonly=TRUE)
par(lty=2, pch=17)
plot(dose, drugA, type="b")
par(opar)

3.1 图例符号和线条

参数 描述 图例
pch 指定绘制点时使用的符号 image
cex 指定符号的大小。cex 是一个数值，表示绘图符号相对于默认大小的缩放倍数。默认大小为 1，1.5 表示放大为默认值的 1.5 倍，0.5 表示缩小为默认值的 50%，等等
lty 指定线条类型 image
lwd 指定线条宽度。lwd 是以默认值的相对大小来表示的（默认值为 1）。例如，lwd=2 将生成一条两倍于默认宽度的线条

3.2 图形颜色

参数 描述
col 默认的绘图颜色。某些函数（如 lines 和 pie）可以接受一个含有颜色值的向量并自动循环使用。例如，如果设定 col=c("red", "blue")并需要绘制三条线，则第一条线将为红色，第二条线为蓝色，第三条线又将为红色
col.axis 坐标轴刻度文字的颜色
col.lab 坐标轴标签（名称）的颜色
col.main 标题颜色
col.sub 副标题颜色
fg 图形的前景色
bg 图形的背景色

创建连续型颜色向量的函数，包括rainbow()、heat.colors()、terrain.colors()、topo.colors()以及cm.colors()

install.packages("RColorBrewer")
library(RColorBrewer) 
n <- 7 
mycolors <- brewer.pal(n, "Set1") 
barplot(rep(1,n), col=mycolors)

# 可选调色板列表
brewer.pal.info
display.brewer.all()

n <- 10 
mycolors <- rainbow(n) 
pie(rep(1, n), labels=mycolors, col=mycolors) 
mygrays <- gray(0:n/n) 
pie(rep(1, n), labels=mygrays, col=mygrays)

3.3 文本属性

参数 描述
cex 表示相对于默认大小缩放倍数的数值。默认大小为 1，1.5 表示放大为默认值的 1.5 倍，0.5 表示缩小为默认值的 50%，等等
cex.axis 坐标轴刻度文字的缩放倍数。类似于 cex
cex.lab 坐标轴标签（名称）的缩放倍数。类似于 cex
cex.main 标题的缩放倍数。类似于 cex
cex.sub 副标题的缩放倍数。类似于 cex
font 整数。用于指定绘图使用的字体样式。1=常规，2=粗体，3=斜体，4=粗斜体，5=符号字体（以 Adobe符号编码表示）
font.axis 坐标轴刻度文字的字体样式
font.lab 坐标轴标签（名称）的字体样式
font.main 标题的字体样式
font.sub 副标题的字体样式
ps 字体磅值（1 磅约为 1/72 英寸）。文本的最终大小为 ps × cex
family 绘制文本时使用的字体族。标准的取值为 serif（衬线）、sans（无衬线）和 mono（等宽）
pin 以英寸表示的图形尺寸（宽和高）
mai 以数值向量表示的边界大小，顺序为“下、左、上、右”，单位为英寸
mar 以数值向量表示的边界大小，顺序为“下、左、上、右”，单位为英分①。默认值为 c(5, 4, 4, 2) + 0.1 
par(font.lab=3, cex.lab=1.5, font.main=4, cex.main=2)

windowsFonts( 
 A=windowsFont("Arial Black"), 
 B=windowsFont("Bookman Old Style"), 
 C=windowsFont("Comic Sans MS") 
)

opar <- par(no.readonly=TRUE) 
par(pin=c(2, 3)) 
par(lwd=2, cex=1.5) 
par(cex.axis=.75, font.axis=3) 
plot(dose, drugA, type="b", pch=19, lty=2, col="red") 
plot(dose, drugB, type="b", pch=23, lty=6, col="blue", bg="green") 
par(opar)

# main-标题，sub-副标题，xlab/ylab-坐标轴标签，xlim/ylim-坐标轴范围
# 高级绘图已经包含了默认的标题和表现，plot()或者par()语句中添加ann=FALSE来移除
plot(dose, drugA, type="b", col="red", lty=2, pch=2, lwd=2, main="Clinical Trials for Drug A", sub="This is hypothetical data", xlab="Dosage", ylab="Drug Response", xlim=c(0, 60), ylim=c(0, 70))

# title()为图形添加标题和坐标轴标签
title(main="main title", sub="subtitle", xlab="x-axis label", ylab="y-axis label")

# 坐标轴axis()，格式为axis(side, at=, labels=, pos=, lty=, col=, las=, tck=, ...)
选项 描述
side 一个整数，表示在图形的哪边绘制坐标轴（1=下，2=左，3=上，4=右）
at 一个数值型向量，表示需要绘制刻度线的位置
labels 一个字符型向量，表示置于刻度线旁边的文字标签（如果为 NULL，则将直接使用 at 中的值）
pos 坐标轴线绘制位置的坐标（即与另一条坐标轴相交位置的值）
lty 线条类型
col 线条和刻度线颜色
las 标签是否平行于（=0）或垂直于（=2）坐标轴
tck 刻度线的长度，以相对于绘图区域大小的分数表示（负值表示在图形外侧，正值表示在图形内侧，0表示禁用刻度，1 表示绘制网格线）；默认值为–0.01 
# 创建自定义坐标轴，应禁用高级绘图函数自动生成的坐标轴
# axes=FALSE将禁用全部坐标轴（包括坐标轴框架线，除非添加了frame.plot=TRUE）
# xaxt="n"和yaxt="n"分别禁用X轴和Y轴（会留下框架先，只除刻度）

# 生成数据
x <- c(1:10)
y <- x
z <- 10/x
opar <- par(no.readonly=TRUE)

# 增加边界大小
par(mar=c(5, 4, 4, 8) + 0.1)
# 绘制x对y的图形
plot(x, y, type="b",
 pch=21, col="red", 
 yaxt="n", lty=3, ann=FALSE) 
# 添加x对1/x的直线
lines(x, z, type="b", pch=22, col="blue", lty=2)

# 绘制自己的坐标轴
axis(2, at=x, labels=x, col.axis="red", las=2)
axis(4, at=z, labels=round(z, digits=2), col.axis="blue", las=2, cex.axis=0.7, tck=-.01)
# 添加标题和文本
mtext("y=1/x", side=4, line=3, cex.lab=1, las=2, col="blue")
title("An Example of Creative Axes",xlab="X values", ylab="Y=X")
par(opar)

# 次要刻度线，nx和ny分别指定了X轴和Y轴的次要刻度线划分的区间个数
library(Hmisc)
minor.tick(nx=n, ny=n, tick.ratio=n)
# X轴的每两条主刻度线之间添加1条次要刻度线，并在Y轴的每两条主刻度线之间添加2条次要刻度线：
minor.tick(nx=2, ny=3, tick.ratio=0.5)

# 参考线 abline(h=yvalues, v=xvalues)
# 在y为1、5、7的位置添加了水平实线
abline(h=c(1,5,7))
# 在x为1、3、5、7、9的位置添加了垂直的蓝色虚线
abline(v=seq(1, 10, 2), lty=2, col="blue")

# 图例 legend(location, title, legend, ...)
选项 描述
location 有许多方式可以指定图例的位置。你可以直接给定图例左上角的 x、y 坐标，也可以执行 locator(1)，然后通过鼠标单击给出图例的位置，还可以使用关键字 bottom、bottomleft、left、topleft、top、topright、right、bottomright 或 center 放置图例。如果你使用了以上某个关键字，那么可以同时使用参数 inset=指定图例向图形内侧移动的大小（以绘图区域大小的分数表示）
title 图例标题的字符串（可选）
legend 图例标签组成的字符型向量
... 其他选项。如果图例标示的是颜色不同的线条，需要指定 col=加上颜色值组成的向量。如果图例标示的是符号不同的点，则需指定 pch=加上符号的代码组成的向量。如果图例标示的是不同的线条宽度或线条类型，请使用 lwd=或 lty=加上宽度值或类型值组成的向量。要为图例创建颜色填充的盒形（常见于条形图、箱线图或饼图），需要使用参数 fill=加上颜色值组成的向量 
dose <- c(20, 30, 40, 45, 60) 
drugA <- c(16, 20, 27, 40, 60) 
drugB <- c(15, 18, 25, 31, 40)
opar <- par(no.readonly=TRUE)
# 增加线条、文本、符号、标签的宽度或大小
par(lwd=2, cex=1.5, font.lab=2)

# 绘制图形
plot(dose, drugA, type="b", 
 pch=15, lty=1, col="red", ylim=c(0, 60),
 main="Drug A vs. Drug B", 
 xlab="Drug Dosage", ylab="Drug Response")
lines(dose, drugB, type="b", pch=17, lty=2, col="blue")
abline(h=c(30), lwd=1.5, lty=2, col="gray")

# 添加次要刻度线
minor.tick(nx=3, ny=3, tick.ratio=0.5)
# 添加图例
legend("topleft", inset=.05, title="Drug Type", c("A","B"), lty=c(1, 2), pch=c(15, 17), col=c("red", "blue"))
par(opar)

# text()可向绘图区域内部添加文本
# text(location, "text to place", pos, ...)
# mtext()向图形的四个边界之一添加文本
# mtext("text to place", side, line=n, ...)
选项 描述
location 文本的位置参数。可为一对 x、y 坐标，也可通过指定 location 为 locator(1)使用鼠标交互式地确定摆放位置
pos 文本相对于位置参数的方位。1=下，2=左，3=上，4=右。如果指定了 pos，就可以同时指定参数 offset=作为偏移量，以相对于单个字符宽度的比例表示
side 指定用来放置文本的边。1=下，2=左，3=上，4=右。你可以指定参数 line=来内移或外移文本，随着值的增加，文本将外移。也可使用 adj=0 将文本向左下对齐，或使用 adj=1 右上对齐 
attach(mtcars)
plot(wt, mpg, main="Mileage vs. Car Weight", xlab="Weight", ylab="Mileage", pch=18, col="blue")
text(wt, mpg, row.names(mtcars), cex=0.6, pos=4, col="red") 
detach(mtcars)

# 显示不同的字体
opar <- par(no.readonly=TRUE) 
par(cex=1.5) 
plot(1:7,1:7,type="n") 
text(3,3,"Example of default text") 
text(4,4,family="mono","Example of mono-spaced text") 
text(5,5,family="serif","Example of serif text") 
par(opar)

# par()函数中使用图形参数mfrow=c(nrows, ncols)来创建按行填充的、行数为nrows、列数为ncols的图形矩阵。
# par(mfrow=c(2,2)),四个图像，2×2
attach(mtcars) 
opar <- par(no.readonly=TRUE) 
par(mfrow=c(2,2))
plot(wt,mpg, main="Scatterplot of wt vs. mpg") 
plot(wt,disp, main="Scatterplot of wt vs. disp") 
hist(wt, main="Histogram of wt") 
boxplot(wt, main="Boxplot of wt") 
par(opar) 
detach(mtcars)

# par(mfrow=c(3,1))，三个图像，3×1
attach(mtcars) 
opar <- par(no.readonly=TRUE) 
par(mfrow=c(3,1)) 
hist(wt) 
hist(mpg) 
hist(disp) 
par(opar) 
detach(mtcars)

# 三个图像，2×2布局，1 1 2 3，第一幅图占横向两个位置
attach(mtcars) 
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE)) 
hist(wt) 
hist(mpg) 
hist(disp) 
detach(mtcars)

# 三个图像，2×2布局，第一行高度：第二行高度=1:2；第二行第一个图像宽度：第二行第二个图像宽度=3:1
attach(mtcars) 
layout(matrix(c(1, 1, 2, 3), 2, 2, byrow = TRUE), 
 widths=c(3, 1), heights=c(1, 2)) 
hist(wt) 
hist(mpg) 
hist(disp) 
detach(mtcars)

# 边界上添加了两幅箱线图的散点图
opar <- par(no.readonly=TRUE) 
par(fig=c(0, 0.8, 0, 0.8)) 
plot(mtcars$wt, mtcars$mpg, 
 xlab="Miles Per Gallon", 
 ylab="Car Weight") 
par(fig=c(0, 0.8, 0.55, 1), new=TRUE) 
boxplot(mtcars$wt, horizontal=TRUE, axes=FALSE)
par(fig=c(0.65, 1, 0, 0.8), new=TRUE) 
boxplot(mtcars$mpg, axes=FALSE)
mtext("Enhanced Scatterplot", side=3, outer=TRUE, line=-3) 
par(opar)

chapter 4 基本数据管理

# 创建leadership数据框
manager <- c(1, 2, 3, 4, 5) 
date <- c("10/24/08", "10/28/08", "10/1/08", "10/12/08", "5/1/09") 
country <- c("US", "US", "UK", "UK", "UK") 
gender <- c("M", "F", "F", "M", "F") 
age <- c(32, 45, 25, 39, 99) 
q1 <- c(5, 3, 3, 3, 2) 
q2 <- c(4, 5, 5, 3, 2) 
q3 <- c(5, 2, 5, 4, 1) 
q4 <- c(5, 5, 5, NA, 2) 
q5 <- c(5, 5, 2, NA, 1) 
leadership <- data.frame(manager, date, country, gender, age, 
 q1, q2, q3, q4, q5, stringsAsFactors=FALSE)

# 将两个新变量整合到原始的数据框中，有三种方式
# 方式1
mydata<-data.frame(x1 = c(2, 2, 6, 4), 
 x2 = c(3, 4, 2, 8))
mydata$sumx <- mydata$x1 + mydata$x2 
mydata$meanx <- (mydata$x1 + mydata$x2)/2
# 方式2
attach(mydata) 
mydata$sumx <- x1 + x2 
mydata$meanx <- (x1 + x2)/2 
detach(mydata)
# 方式3 推荐
mydata <- transform(mydata, sumx = x1 + x2, meanx = (x1 + x2)/2)

# 将99岁的年龄值重编码为缺失值
leadership$age[leadership$age == 99] <- NA
# variable[condition] <- expression将仅在condition的值为TRUE时执行赋值
leadership$agecat[leadership$age > 75] <- "Elder" 
leadership$agecat[leadership$age >= 55 & 
 leadership$age <= 75] <- "Middle Aged" 
leadership$agecat[leadership$age < 55] <- "Young"
# or
leadership <- within(leadership,{ 
 agecat <- NA 
 agecat[age > 75] <- "Elder" 
 agecat[age >= 55 & age <= 75] <- "Middle Aged" 
 agecat[age < 55] <- "Young" })
 
# 变量的重命名fix(leadership)，交互式修改
# names()函数来重命名
names(leadership)[2] <- "testData"
names(leadership)[6:10] <- c("item1", "item2", "item3", "item4", "item5")

# plyr包有个rename()函数，可用于修改变量名
# rename()格式rename(dataframe, c(oldname="newname", oldname="newname",...))
install.packages("plyr")
library(plyr)
leadership <- rename(leadership, c(manager="managerID", date="testDate"))

# 判断是否为缺失值，is.na(y)，NA->not available
# 缺失值不可比较；缺失值不是无限或不可能出现的数值，需要用到is.infinite() or is.nan()
is.na(leadership[,6:10])

# na.rm=TRUE，排除掉缺失值并进行分析
x <- c(1, 2, NA, 3)
y <- sum(x, na.rm=TRUE)
# y = 6

# na.omit()移除所有含有缺失值的观测
newdata <- na.omit(leadership)
newdata

日期值

符号 含义 示例
%d 数字表示的日期（0-31） 01~31
%a 缩写的星期名 Mon
%A 非缩写的星期名 Monday
%m 月份（00~12） 00~12
%b 缩写的月份 Jan
%B 非缩写的月份 January
%y 两位数的年份 23
%Y 四位数的年份 2023 
# Sys.Data()可以返回当天的日期
# data()返回当天的日期和时间
# format(x, format="output_format")来输出指定格式的日期值，或提取部分数据
today <- Sys.Date()
format(today, format="%B %d %Y")
format(today, format="%A")

# R存储的日期可进行算术计算
startdate <- as.Date("2004-02-13")
enddate <- as.Date("2011-01-22")
days <- enddate - startdate
days

today <- Sys.Date()
dob <- as.Date("1956-10-12")
difftime(today, dob, units="weeks")

# 将日期转换成字符型变量as.character()
strDates <- as.character(dates)

类型转换

判断 转换
is.numeric() as.numeric()
is.character() as.character()
is.vector() as.vector()
is.matrix() as.matrix()
is.data.frame() as.data.frame()
is.factor() as.factor()
is.logical() as.logical() 
# order()函数对一个数据框进行排序
# 按经理人的性别和年龄降序排序
attach(leadership)
newdata <-leadership[order(gender, -age),]
detach(leadership)

# 数据集的合并
# 1、添加列，横向合并两个数据框(集)，使用merge()，即inner join，eg. total <- merge(dataframeA, dataframeB, by="ID")
# or total <- merge(dataframeA, dataframeB, by=c("ID","Country"))
# 使用cbind()进行横向合并，前提是拥有相同的行数，eg. total <- cbind(A, B)
# 2、添加行，total <- rbind(A, B)，必须拥有相同的变量。或删除或追加为NA

# 选取子集，即切片，使用过dataframe[row indices, column indices]
newdata <- leadership[, c(6:10)]

myvars <- c("q1", "q2", "q3", "q4", "q5")
newdata <-leadership[myvars]
# 等价于
myvars <- paste("q", 1:5, sep="")
newdata <-leadership[myvars]

# 丢弃变量
# 知道序号，q3 <=> 8，newdata <- leadership[c(-8,-9)]
myvars <- names(leadership) %in% c("q3", "q4")
newdata <- leadership[!myvars]

# 使用subset()函数进行数据切片并进行处理
newdata <- subset(leadership, age >= 35 | age < 24, select=c(q1, q2, q3, q4))

# 随机抽样
# sampling包抽取和校正调查样本
# survery包分析复杂调查数据
# sample()有放回或无放回抽取大小为n的随机样本
mysample <- leadership[sample(1:nrow(leadership), 3, replace=FALSE),]

# SQL相关
install.packages("sqldf")
library(sqldf)
# mpg对车型进行升序排列，并保留原始行名1
newdf <- sqldf("select * from mtcars where carb=1 order by mpg", row.names=TRUE)
# 输出四缸和六缸车型每一gear水平的mpg和disp的平均值
sqldf("select avg(mpg) as avg_mpg, avg(disp) as avg_disp, gear 
 from mtcars where cyl in (4, 6) group by gear")

chapter 5 高级数据管理

数学函数

函数 描述
abs(x) 绝对值
sqrt(x) 平方根，与x^(0.5)等价
ceiling(x) 不小于x的最小整数
floor(x) 不大于x的最大整数
trunc(x) 向0 的方向截取的 x 中的整数部分
round(x,digits=n) 将 x 舍入为指定位的小数
signif(x,digits=n) 将 x 舍入为指定的有效数字位数
cos(x)、sin(x)、tan(x) 余弦、正弦和正切
acos(x)、asin(x)、atant(x) 反余弦、反正弦和反正切
cosh(x)、sinh(x)、tanh(x) 双曲余弦、双曲正弦和双曲正切
acosh(x)、asinh(x)、atanh(x) 反双曲余弦、反双曲正弦和反双曲正切
log(x,base=n) 对x取以n为底的对数；log(x)自然对数；log10(x)为常用对数
exp(x) 指数函数

统计函数

函数 描述
mean(x) 平均数
median(x) 中位数
sd(x) 标准差
var(x) 方差
mad(x) 绝对中位差（median absolute deviation）
quantile(x,probs) 求分位数。其中 x 为待求分位数的数值型向量，probs 为一个由[0,1]之间的概率值组成的数值向量
range(x) 求值域
sum(x) 求和
diff(x, lag=n) 滞后差分，lag 用以指定滞后几项。默认的 lag 值为 1
min(x) 最小值
max(x) 最大值
scale(x,center=TRUE, scale=TRUE) 为数据对象 x 按列进行中心化(center=TRUE)或标准化(center=TRUE,scale=TRUE)； 
# 函数scale()对矩阵或数据框的指定列进行均值为0、标准差为1的标准化
newdata <- scale(mydata)
# 每一列进行任意均值和标准差的标准化
newdata <- scale(mydata)*SD + M
# 对指定列进行标准化
newdata <- transform(mydata, myvar = scale(myvar)*10+50)

概率函数
[dpqr]distribution_abbreviation()

  • d = 密度函数（density）
  • p = 分布函数（distribution function）
  • q = 分位数函数（quantile function）
  • r = 生成随机数（随机偏差）
分布名称 缩写 分部名称 缩写
Beta 分布 beta Logistic 分布 logis
二项分布 binom 多项分布 multinom
柯西分布 cauchy 负二项分布 nbinom
（非中心）卡方分布 chisq 正态分布 norm
指数分布 exp 泊松分布 pois
F分布 f Wilcoxon 符号秩分布 signrank
Gamma 分布 gamma t 分布 t
几何分布 geom 均匀分布 unif
超几何分布 geom Weibull 分布 weibull
对数正态分布 lnorm Wilcoxon 秩和分布 wilcox 
# 在区间[–3, 3]上绘制标准正态曲线
x <- pretty(c(-3,3), 30) 
y <- dnorm(x) 
plot(x, y, 
 type = "l", 
 xlab = "Normal Deviate", 
 ylab = "Density", 
 yaxs = "i"
)

# 密度函数（dnorm）、分布函数（pnorm）、分位数函数（qnorm）和随机数生成函数（rnorm）
# 位于 z=1.96 左侧的标准正态曲线下方面积
pnorm(1.96)
# 均值为 500，标准差为 100 的正态分布的 0.9 分位点值
qnorm(.9, mean=500, sd=100)
# 生成 50 个均值为 50，标准差为 10 的正态随机数
rnorm(50, mean=50, sd=10)

# 通过set.seed()显式指定种子，让结果重现（reproducible）
# 生成自给定均值向量和协方差阵的多元正态数据，使用MASS包中的mvrnorm()函数，mvrnorm(n, mean, sigma)
# n是想要的样本大小，mean为均值向量，而sigma是方差协方差矩阵（或相关矩阵）
library(MASS)
options(digits=3)
# 设定随机数种子
set.seed(1234)
# 指定均值向量、协方差阵
mean <- c(230.7, 146.7, 3.6)
sigma <- matrix(c(15360.8, 6721.2, -47.1, 
 6721.2, 4700.9, -16.5, 
 -47.1, -16.5, 0.3), nrow=3, ncol=3)
# 生成数据
mydata <- mvrnorm(500, mean, sigma)
mydata <- as.data.frame(mydata)
names(mydata) <- c("y","x1","x2")

# 查看结果
dim(mydata)
head(mydata, n=10)

字符处理函数

函数 描述
nchar(x) 计算 x 中的字符数量
substr(x,start,stop) 提取或替换一个字符向量中的子串
grep(pattern,x,ignore.case=FALSE, fixed=FALSE) 在 x 中搜索某种模式。若 fixed=FALSE，则 pattern 为一个正则表达式。若fixed=TRUE，则 pattern 为一个文本字符串。返回值为匹配的下标
sub(pattern, replacement, x, ignore.case=FALSE, fixed=FALSE) 在 x 中搜索 pattern，并以文本 replacement 将其替换。若 fixed=FALSE，则pattern 为一个正则表达式。若 fixed=TRUE，则 pattern 为一个文本字符串。
strsplit(x, split, fixed=FALSE) 在 split 处分割字符向量 x 中的元素。若 fixed=FALSE，则 pattern 为一个正则表达式。若 fixed=TRUE，则 pattern 为一个文本字符串
paste(…, sep="") 连接字符串，分隔符为 sep；paste("x", 1:3,sep="")返回值为 c("x1", "x2", "x3")
toupper(x) 大写转换
tolower(x) 小写转换
length(x) 对象 x 的长度
seq(from, to, by) 生成一个序列；indices <- seq(1,10,2)，indices 的值为 c(1, 3, 5, 7, 9)
rep(x, n) 将x 重复 n 次
cut(x, n) 将连续型变量 x 分割为有着 n 个水平的因子，使用选项 ordered_result = TRUE 以创建一个有序型因子
pretty(x, n) 创建美观的分割点。通过选取 n+1 个等间距的取整值，将一个连续型变量 x分割为 n 个区间。
cat(... , file ="myfile", append =FALSE) 连接...中的对象，并将其输出到屏幕上或文件中（如果声明了一个的话）

\n表示新行，\t为制表符，'为单引号，\b为退格，等等。（键入?Quotes以了解更多。）

# apply()函数，可将一个任意函数“应用”到矩阵、数组、数据框的任何维度上。apply()函数的使用格式为：apply(x, MARGIN, FUN, ...)
mydata <- matrix(rnorm(30), nrow=6)

# 计算每行均值
apply(mydata, 1, mean)
# 计算每列均值
apply(mydata, 2, mean)
# 计算每列的截尾均值，取中间60%的数据，首尾各20%的数据舍去
apply(mydata, 2, mean, trim=0.2)

# 数据处理例子
# 限制小数点后两位，更易阅读
options(digits=2)
# 输入数据
Student <- c("John Davis", "Angela Williams", "Bullwinkle Moose", 
 "David Jones", "Janice Markhammer", "Cheryl Cushing", 
 "Reuven Ytzrhak", "Greg Knox", "Joel England", 
 "Mary Rayburn")
Math <- c(502, 600, 412, 358, 495, 512, 410, 625, 573, 522)
Science <- c(95, 99, 80, 82, 75, 85, 80, 95, 89, 86)
English <- c(25, 22, 18, 15, 20, 28, 15, 30, 27, 18)
roster <- data.frame(Student, Math, Science, English, stringsAsFactors=FALSE)
# 计算综合得分，先进行标准化
z <- scale(roster[,2:4])
# 添加均值计算综合得分
score <- apply(z, 1, mean)
# 使用cbind进行综合得分列添加
roster <- cbind(roster, score)
# 使用综合分数的百分位进行分级判定
y <- quantile(score, c(.8,.6,.4,.2))
# 使用逻辑运算符进行分级判定
roster$grade[score >= y[1]] <- "A"
roster$grade[score < y[1] & score >= y[2]] <- "B"
roster$grade[score < y[2] & score >= y[3]] <- "C"
roster$grade[score < y[3] & score >= y[4]] <- "D"
roster$grade[score < y[4]] <- "F"
# 使用空格拆分姓名
name <- strsplit((roster$Student), " ")
# 使用sapply()提取姓和名
Lastname <- sapply(name, "[", 2) 
Firstname <- sapply(name, "[", 1)
# 加入分开的姓和名，并舍弃之前的Student列
roster <- cbind(Firstname,Lastname, roster[,-1])
# 对新数据进行排序
roster <- roster[order(Lastname,Firstname),]

控制流

  • 语句statement
  • 条件cond
  • 表达式expr
  • 序列seq
# 1 循环
# for (var in seq) statement
for (i in 1:10) print("Hello")
# while (cond) statement
i <- 10
while (i > 0) {print("Hello"); i <- i - 1}
# 2 条件执行
# if (cond) statement
# if (cond) statement1 else statement2
if (is.character(grade)) grade <- as.factor(grade) 
if (!is.factor(grade)) grade <- as.factor(grade) else print("Grade already is a factor")
# ifelse(cond, statement1, statement2)
ifelse(score > 0.5, print("Passed"), print("Failed"))
outcome <- ifelse (score > 0.5, "Passed", "Failed")
# switch(expr, ...)
feelings <- c("sad", "afraid")
for (i in feelings)
print(
switch(i,
happy  = "I am glad you are happy",
afraid = "There is nothing to fear",
sad    = "Cheer up",
angry  = "Calm down now"
)
)

# 自定义函数
mystats <- function(x, parametric=TRUE, print=FALSE) { 
 if (parametric) { 
 center <- mean(x); spread <- sd(x) 
 } else { 
 center <- median(x); spread <- mad(x) 
 } 
 if (print & parametric) { 
 cat("Mean=", center, "\n", "SD=", spread, "\n") 
 } else if (print & !parametric) { 
 cat("Median=", center, "\n", "MAD=", spread, "\n") 
 } 
 result <- list(center=center, spread=spread) 
 return(result) 
}
set.seed(1234) 
x <- rnorm(500)
y <- mystats(x)
y <- mystats(x, parametric=FALSE, print=TRUE)

mydate <- function(type="long") { 
 switch(type, 
 long = format(Sys.time(), "%A %B %d %Y"), 
 short = format(Sys.time(), "%m-%d-%y"), 
 cat(type, "is not a recognized type\n") 
 ) 
}
mydate("long")
mydate("short")
mydate()
mydate("medium")
# 函数cat()仅会在输入不匹配时执行
# 函数warning()来生成一条错误提示信息
# 函数用message()来生成一条诊断信息
# 函数stop()停止当前表达式的执行并提示错误
# 可参考Venables & Ripley（2000）以及Chambers（2008）。

# 转置
cars <- mtcars[1:5, 1:4]
cars
t(cars)
# 整合aggregate(x, by, FUN)，x是待折叠的数据对象，by是一个变量名组成的列表，而FUN则是用来计算描述性统计量的标量函数
aggdata <-aggregate(mtcars, by=list(cyl,gear), FUN=mean, na.rm=TRUE)
# reshape2包是一套重构和整合数据集的绝妙的万能工具
install.packages("reshape2")
library(reshape2)
ID <- c(1,1,2,2)
Time <- c(1,2,1,2)
X1 <- c(5,3,6,2)
X2 <- c(6,5,1,4)
mydata <- data.frame(ID, Time, X1, X2)
# melt()融合，每个测量变量独占一行，指X1和X2
md <- melt(mydata,id=c("ID","Time"))
# dcast()重铸newdata <- dcast(md, formula, fun.aggregate)

image

chapter 6 基本图形

# 条形图barplot(height)
install.packages("vcd")
library(vcd)
counts <- table(Arthritis$Improved)
# 简单条形图
barplot(counts, 
 main="Simple Bar Plot", 
 xlab="Improvement", ylab="Frequency")
# 水平条形图
barplot(counts, 
 main="Horizontal Bar Plot", 
 xlab="Frequency", ylab="Improvement",
 horiz=TRUE)
 
counts <- table(Arthritis$Improved, Arthritis$Treatment)
# 堆砌条形图
barplot(counts, 
 main="Stacked Bar Plot", 
 xlab="Treatment", ylab="Frequency", 
 col=c("red", "yellow","green"), 
 legend=rownames(counts))
# 分组条形图
barplot(counts, 
 main="Grouped Bar Plot", 
 xlab="Treatment", ylab="Frequency", 
 col=c("red", "yellow", "green"), 
 legend=rownames(counts), beside=TRUE)
 
# 均值条形图
states <- data.frame(state.region, state.x77)
means <- aggregate(states$Illiteracy, by=list(state.region), FUN=mean)
# 均值排序
means <- means[order(means$x),]
barplot(means$x, names.arg=means$Group.1)

# 条形图微调，搭配标签
# 增加y边界的大小
par(mar=c(5,8,4,2))
# 旋转条形的标签
par(las=2)
counts <- table(Arthritis$Improved)
# 缩小字体大小，修改标签文本
barplot(counts, 
 main="Treatment Outcome", 
 horiz=TRUE, 
 cex.names=0.8,
 names.arg=c("No Improvement", "Some Improvement", 
 "Marked Improvement"))

# 棘状图spinogram，对堆砌条形图进行了重缩放，类似于标准化
library(vcd) 
attach(Arthritis) 
counts <- table(Treatment, Improved) 
spine(counts, main="Spinogram Example") 
detach(Arthritis)

# 饼图pie(x, labels)----不建议，看个人取舍
install.packages("plotrix")
par(mfrow=c(2, 2))
slices <- c(10, 12,4, 16, 8) 
lbls <- c("US", "UK", "Australia", "Germany", "France") 
pie(slices, labels = lbls, 
 main="Simple Pie Chart") 
pct <- round(slices/sum(slices)*100) 
lbls2 <- paste(lbls, " ", pct, "%", sep="") 
pie(slices, labels=lbls2, col=rainbow(length(lbls2)),
 main="Pie Chart with Percentages") 
library(plotrix) 
pie3D(slices, labels=lbls,explode=0.1, 
 main="3D Pie Chart ") 
mytable <- table(state.region)
lbls3 <- paste(names(mytable), "\n", mytable, sep="") 
pie(mytable, labels = lbls3, 
 main="Pie Chart from a Table\n (with sample sizes)")
 
# 扇形图
library(plotrix) 
slices <- c(10, 12,4, 16, 8) 
lbls <- c("US", "UK", "Australia", "Germany", "France") 
fan.plot(slices, labels = lbls, main="Fan Plot")

# 直方图hist(x)
par(mfrow=c(2,2)) 
# 简单直方图
hist(mtcars$mpg)
# 指定组数和颜色
hist(mtcars$mpg, 
 breaks=12, 
 col="red", 
 xlab="Miles Per Gallon", 
 main="Colored histogram with 12 bins")
# 添加轴须图
hist(mtcars$mpg, 
 freq=FALSE, 
 breaks=12, 
 col="red", 
 xlab="Miles Per Gallon", 
 main="Histogram, rug plot, density curve")
rug(jitter(mtcars$mpg)) 
lines(density(mtcars$mpg), col="blue", lwd=2)

# 添加正态密度曲线和外框
x <- mtcars$mpg 
h<-hist(x, 
 breaks=12, 
 col="red", 
 xlab="Miles Per Gallon", 
 main="Histogram with normal curve and box") 
xfit<-seq(min(x), max(x), length=40) 
yfit<-dnorm(xfit, mean=mean(x), sd=sd(x)) 
yfit <- yfit*diff(h$mids[1:2])*length(x)
lines(xfit, yfit, col="blue", lwd=2) 
box()

# 核密度plot(density(x))
par(mfrow=c(2,1)) 
d <- density(mtcars$mpg) 
# 使用默认配置生成核密度图
plot(d)
d <- density(mtcars$mpg)
# 添加一个标题
plot(d, main="Kernel Density of Miles Per Gallon")
# 将曲线修改为蓝色，并使用实心红色填充曲线下方的区域
polygon(d, col="red", border="blue")
# 添加棕色的轴须图
rug(mtcars$mpg, col="brown")

# 使用sm包中的sm.density.compare()函数可向图形叠加两组或更多的核密度图。sm.density.compare(x, factor)
# 可比较的核密度图
library(sm) 
attach(mtcars) 
# 创建分组因子
cyl.f <- factor(cyl, levels= c(4,6,8), 
 labels = c("4 cylinder", "6 cylinder",
 "8 cylinder"))
# 绘制密度图
sm.density.compare(mpg, cyl, xlab="Miles Per Gallon") 
title(main="MPG Distribution by Car Cylinders")
# 通过鼠标点击添加图例
colfill<-c(2:(1+length(levels(cyl.f)))) 
legend(locator(1), levels(cyl.f), fill=colfill)
detach(mtcars)

# 箱线图（盒须图）boxplot(mtcars$mpg, main="Box plot", ylab="Miles per Gallon")

image

# 箱线图展示单个变量或分组变量boxplot(formula, data=dataframe)
boxplot(mpg ~ cyl, data=mtcars, 
 main="Car Mileage Data", 
 xlab="Number of Cylinders", 
 ylab="Miles Per Gallon")
 
# 含凹槽的箱线图
boxplot(mpg ~ cyl, data=mtcars, 
 notch=TRUE, 
 varwidth=TRUE, 
 col="red", 
 main="Car Mileage Data", 
 xlab="Number of Cylinders", 
 ylab="Miles Per Gallon")
 
# 两个交叉因子的箱线图
# 创建汽缸数量的因子
mtcars$cyl.f <- factor(mtcars$cyl, 
 levels=c(4,6,8), 
 labels=c("4","6","8"))
# 创建变速箱类型的因子
mtcars$am.f <- factor(mtcars$am, 
 levels=c(0,1), 
 labels=c("auto", "standard"))
# 生成箱线图
boxplot(mpg ~ am.f *cyl.f, 
 data=mtcars, 
 varwidth=TRUE, 
 col=c("gold","darkgreen"), 
 main="MPG Distribution by Auto Type", 
 xlab="Auto Type", ylab="Miles Per Gallon")

# 小提琴图
library(zoo)
library(vioplot)
x1 <- mtcars$mpg[mtcars$cyl==4] 
x2 <- mtcars$mpg[mtcars$cyl==6] 
x3 <- mtcars$mpg[mtcars$cyl==8] 
vioplot(x1, x2, x3, 
 names=c("4 cyl", "6 cyl", "8 cyl"), 
 col="gold") 
title("Violin Plots of Miles Per Gallon", ylab="Miles Per Gallon", 
 xlab="Number of Cylinders")

# 点图dotchart(x, labels=)
dotchart(mtcars$mpg, labels=row.names(mtcars), cex=.7, 
 main="Gas Mileage for Car Models", 
 xlab="Miles Per Gallon")

# 分组->排序->着色点图
# 每加仑汽油书进行排序后保存为数据框x
x <- mtcars[order(mtcars$mpg),]
# 转换成因子
x$cyl <- factor(x$cyl)
# 添加字符型向量分别进行着色
x$color[x$cyl==4] <- "red"
x$color[x$cyl==6] <- "blue"
x$color[x$cyl==8] <- "darkgreen"
# 进行分组，添加颜色等
dotchart(x$mpg, 
 labels = row.names(x),
 cex=.7, 
 groups = x$cyl, 
 gcolor = "black", 
 color = x$color,
 pch=19, 
 main = "Gas Mileage for Car Models\ngrouped by cylinder", 
 xlab = "Miles Per Gallon")

chapter 7 基本统计分析

可参考文献McCall（2000）和Kirk（2007）

myvars <- c("mpg", "hp", "wt")
# 显示前6行数据
head(mtcars[myvars])
# 显示描述统计量
summary(mtcars[myvars])
# 通过sapply()计算描述性统计量，包括偏度和峰度
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, stdev=s, skew=skew, kurtosis=kurt)) 
 }
sapply(mtcars[myvars], mystats)
# 单纯地忽略缺失值，那么应当使用sapply(mtcars[myvars], mystats, na.omit=TRUE)

# 其他包计算描述统计量，包括Hmisc、pastecs和psych
# Hmisc包中的describe()函数可返回变量和观测的数量、缺失值和唯一值的数目、平均值、分位数，以及五个最大的值和五个最小的值
library(Hmisc)
describe(mtcars[myvars])
# pastecs包中有一个名为stat.desc()的函数，stat.desc(x, basic=TRUE, desc=TRUE, norm=FALSE, p=0.95)
# x是一个数据框或时间序列
# basic=TRUE（默认值），则计算其中所有值、空值、缺失值的数量，以及最小值、最大值、值域，还有总和
# desc=TRUE（默认值），则计算中位数、平均数、平均数的标准误、平均数置信度为95%的置信区间、方差、标准差以及变异系数
# norm=TRUE（非默认值），则返回正态分布统计量，包括偏度和峰度（以及它们的统计显著程度）和Shapiro-Wilk正态检验结果
# p值来计算平均数的置信区间

# 通过pastecs包中的stat.desc()函数计算描述性统计量
install.packages("pastecs")
library(pastecs)
stat.desc(mtcars[myvars])
# 通过psych包中的describe()计算描述性统计量
install.packages("psych")
library(psych)
describe(mtcars[myvars])

# 使用aggregate()分组获取描述性统计量
aggregate(mtcars[myvars], by=list(am=mtcars$am), mean)
aggregate(mtcars[myvars], by=list(am=mtcars$am), sd)
# aggregate()仅允许在每次调用中使用平均数、标准差这样的单返回值函数。它无法一次返回若干个统计量。要完成这项任务，可以使用by()函数。格式为：by(data, INDICES, FUN)
dstats <- function(x)sapply(x, mystats)
by(mtcars[myvars], mtcars$am, dstats)

# doBy包和psych包也提供了分组计算描述性统计量的函数
install.packages("doBy")
# 如遇版本过低，则在文件夹里搜索相关包，手动删除后再安装，成功（需要退出R软件）
library(doBy)
summaryBy(mpg+hp+wt~am, data=mtcars, FUN=mystats)
# 使用psych包中的describeBy()分组计算概述统计量
myvars <- c("mpg", "hp", "wt")
describeBy(mtcars[myvars], list(am=mtcars$am))

频数

函数 描述
table(var1, var2, ..., varN) 使用 N 个类别型变量（因子）创建一个 N 维列联表
xtabs(formula, data) 根据一个公式和一个矩阵或数据框创建一个 N 维列联表
prop.table(table, margins) 依 margins 定义的边际列表将表中条目表示为分数形式
margin.table(table, margins) 依 margins 定义的边际列表计算表中条目的和
addmargins(table, margins) 将概述边 margins（默认是求和结果）放入表中
ftable(table) 创建一个紧凑的“平铺”式列联表 
# 一维列联表
mytable <- with(Arthritis, table(Improved))
# 用prop.table()将这些频数转化为比例值；prop.table()*100转换成百分比
prop.table(mytable)

# 二维列联表：mytable <- table(A, B)；mytable <- xtabs(~ A + B, data=mydata)；用gmodels包中的CrossTable()函数
mytable <- xtabs(~ Treatment+Improved, data=Arthritis)
# 用margin.table()和prop.table()函数分别生成边际频数和比例。
# 下标1指代table()语句中的第一个变量。
margin.table(mytable, 1)
prop.table(mytable, 1)
# 用addmargins()函数为这些表格添加边际和
addmargins(mytable)
# CrossTable()函数仿照SAS中PROC FREQ或SPSS中CROSSTABS的形式生成二维列联表
library(gmodels) 
CrossTable(Arthritis$Treatment, Arthritis$Improved)

# 多维列联表：使用table() 和 xtabs() 都可以基于三个或更多的类别型变量生成多维列联表。
mytable <- xtabs(~ Treatment+Sex+Improved, data=Arthritis)
# ftable()函数可以以一种紧凑而吸引人的方式输出多维列联表
ftable(mytable)
# 边际频数
margin.table(mytable, 1)
margin.table(mytable, c(1, 3))
ftable(prop.table(mytable, c(1, 2)))
ftable(addmargins(prop.table(mytable, c(1, 2)), 3))

# 独立性检验
# 使用chisq.test()函数对二维表的行变量和列变量进行卡方独立性检验
library(vcd)
mytable <- xtabs(~Treatment+Improved, data=Arthritis)
chisq.test(mytable)
# p-value = 0.001463 治疗情况和改善情况不独立
mytable <- xtabs(~Improved+Sex, data=Arthritis) 
chisq.test(mytable)
# p-value = 0.08889 性别和改善情况独立

# fisher.test()函数进行Fisher精确检验。
# Fisher精确检验的原假设是：边界固定的列联表中行和列是相互独立的。其调用格式为fisher.test(mytable)，其中的mytable是一个二维列联表。
# fisher.test()函数可以在任意行列数大于等于2的二维列联表上使用，但不能用于2×2的列联表。
mytable <- xtabs(~Treatment+Improved, data=Arthritis)
fisher.test(mytable)

# mantelhaen.test()函数可用来进行Cochran-Mantel-Haenszel卡方检验。
# 原假设是，两个名义变量在第三个变量的每一层中都是条件独立的。下列代码可以检验治疗情况和改善情况在性别的每一水平下是否独立。此检验假设不存在三阶交互作用（治疗情况×改善情况×性别）。
mytable <- xtabs(~Treatment+Improved+Sex, data=Arthritis)
mantelhaen.test(mytable)

# 相关性度量
# vcd包中的assocstats()函数可以用来计算二维列联表的phi系数、列联系数和Cramer’s V系数。
library(vcd)
mytable <- xtabs(~Treatment+Improved, data=Arthritis)
assocstats(mytable)
# vcd包也提供了一个kappa()函数，可以计算混淆矩阵的Cohen’s kappa值以及加权的kappa值。

# 相关：相关系数的符号（±）表明关系的方向（正相关或负相关），其值的大小表示关系的强弱程度（完全不相关时为0，完全相关时为1）

参考资料

[1] John Cook - R programming for those coming from other languages
[2] Google’s R Style Guide

posted @ 2023-03-13 08:55  CestLaVie0019  阅读(189)  评论(0)    收藏  举报