R 语言实战-Part 5-2笔记
R 语言实战（第二版）
----------第21章创建包--------------------------
#包是一套函数、文档和数据的合集，以一种标准的格式保存
#1.测试npar包。进行非参组间比较
pkg <- "npar_1.0.tar.gz"
loc <- "http://www.statmethods.net/RiA"
url <- paste(loc,pkg,sep = "/")
download.file(url,pkg)
install.packages(pkg,repos = NULL,type = "source")
#假设检验
library(npar)
help("life")
head(life)
hist(life$hlef,col = "grey",breaks = 10) #检查正态性：因变量负偏
library(ggplot2)
ggplot(life,aes(region,hlef))+
geom_point(size=3,color="darkgrey")+theme_bw()#检查方差齐性：方差不同
#比较分组
results <- oneway(hlef~region,life) #计算多重比较时，要考虑到alpha膨胀可能性，p.ajust来校正
summary(results)#包括Kruskal-Wallis检验组间差异，总体统计量，成对分组比较(Wilcoxon秩和检验)
plot(results,col="lightblue")
#2.开发npar包
##1）计算统计量函数：oneway.R
#开头注释#'，会被roxygen2包生成包文档
#' @title 非参组间比较
#'
#' @description
#' \code{oneway} 计算非参组间比较，包括综合检验和事后成对组间比较
#'
#' @details
#' 这个函数计算了一个综合Kruskal-Wallis检验，用于检验组别是否相等，接着使用Wilcoxon秩和检验来进行成对比较
#'
#' @param formula
#' @param data
#' @param exact逻辑型变量
#' @param sort逻辑型变量
#' @param 用于调整多重比较的p值的方法,默认"holm"
#' @export
#' @return
#' \item{CALL}
#' \item{data}
#' \item{sumstats}
#' \item{kw}
#' \item{method}
#' \item{wmc}
#' \item{vnames}
#' @author pjx<pjx@bgi.com>
#' @examples
#' results <- oneway(helf~region,life)
#' summary(results)
#' plot(results,col="lightblue)
oneway <- function(formula,data,exact=FALSE,sort=TRUE,
method=c("holm","hochberg","hommel","bonferroni","BH","BY","fdr","none")){
#检查参数
if(missing(formula) || class(formula) != "formula" || length(all.vars(formula)) != 2)
stop("'formula' is missing or incorrect")
method <- match.arg(method)
#设定数据
df <- model.frame(formula,data) #创建一个数据框，第一列为因变量，第二列为分组变量
y <- df[[1]]
g <- as.factor(df[[2]])
vnames <- names(df)
#重新排序
if(sort) g <- reorder(g,y,FUN=median) #对分组变量g按照因变量y的中位数排序
groups <- levels(g)
k <- nlevels(g) #组的数目
#总体统计量
getstats <- function(x)(c(N=length(x),Median=median(x),MAD=mad(x)))
sumstats <- t(aggregate(y,by=list(g),FUN=getstats)[2])
rownames(sumstats) <- c("n","median","mad")
colnames(sumstats) <- groups
#统计检验
kw <- kruskal.test(formula,data)
wmc <- NULL
for(i in 1:(k-1)){
for (j in (i+1):k) {
y1 <- y[g==groups[i]]
y2 <- y[g==groups[j]]
test <- wilcox.test(y1,y2,exact = exact)
r <- data.frame(Group.1=group[i],Group.2=group[j],
W=test$statistic[[1]],p=test$p.value)
wmc <- rbind(wmc,r)
}
}
wmc$p <- p.adjust(wmc$p,method = method)
#返回结果
data <- data.frame(y,g)
names(data) <- vnames
results <- list(CALL=match.call(), #函数调用
data=data,sumstats=sumstats,kw=kw,
method=method,wmc=wmc,vnames=vnames)
class(results) <- c("oneway","list") #设置这个列表的类
return(results)
}
#尽管以上结果列表包含所有信息，但一般不会直接获取单个分量的信息。因此还需创建泛型函数如print/summary/plot等来表达更为具体和有意义的方法
##2）创建函数：print.R
#' @title
#'
#' @description
#' \code{print.oneway}
#'
#' @details
#'
#' @param x
#' @param ...
#' @method print oneway
#' @export
#' @return
#' @author
#' @example
#' results <- oneway(hlef~region,life)
#' print(resluts)
print.oneway <- function(x,...){
#检查输入
if(!inherits(x,"oneway")) #inherits函数确保有这个类
stop("object must be of class 'oneway'")
#打印头部
cat("data:",x$vnames[1],"by",x$vnames[2],"\n\n")
cat("multiple comparisons (Wilcoxon rank sum tests)\n")
cat(paste("probability adjustment = ", x$method,"\n",sep = ""))
#打印表格
print(x$wmc,...)
}
##3)汇总函数：summary.R
summary(results)
#' @title
#'
#' @description
#' \code{summary.oneway}
#' nonparametric analysis.
#'
#' @details
#'
#' @param object
#' @param ...
#' @method summary oneway
#' @export
#' @return
#' @author
#' @example
#' results <- oneway(hlef~region,life)
#' summary(results)
summrary.oneway <- function(object,...){
if(!inherits(object,"oneway"))
stop("object must be of class 'oneway'")
if(!exists("digits")) digits <- 4L
kw <- object$kw
wmc <- object$wmc
cat("data:",object$vnames[1],"on",object$vnames[2],"\n\n")
#Kruskal-Wallis检验
cat("omnibus test\n")
cat(paste("Kruskal-Wallis chi-squared= ",
round(kw$statistic,4),
",df=",round(kw$parameter,3),
",p-value=",format.pval(kw$p.value,digits = digits),
"\n\n",sep = ""))
#描述性统计量
cat("Descriptive statistics\n")
print(object$sumstats,...)
#表格标记
wmc$stars <- " "
wmc$stars[wmc$p<0.1] <- "."
wmc$stars[wmc$p<0.05] <- "*"
wmc$stars[wmc$p<0.01] <- "**"
wmc$stars[wmc$p<0.001] <- "***"
names(wmc)[which(names(wmc)="stars")] <- " "
#成分分组比较
cat("\nmultiple comparisons (Wilcoxon rank sum tests)\n")
cat(paste0("probability ajustment=",object$method,"\n"))
print(wmc,...)
cat("---\nsignif.codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n")
}
##4)绘图函数：plot.R
#注释同上，略
plot.oneway <- function(x,...){
if(!inherits(x,"oneway"))
stop("object must be of class 'oneway'")
#生成箱线图
data <- x$data
y <- data[,1]
g <- data[,2]
stats <- x$sumstats
lbl <- paste("md=",stats[2,],"\nn=",stats[1,],sep = "")
opar <- par(no.readonly = T)
boxplot(y~g,...)
#为图添加标记
abline(h=median(y),lty=2,col="darkgrey")
axis(3,at=1:length(lbl),labels = lbl,cex.axis=0.9)
par(opar)
}
##5）添加demo数据到包
region <- c()
state <- c()
hlem <- c()
hlef <- c()
life <- data.frame(region,state,hlem,hlef)
save(life,file = "life.rda")
#添加数据也需要一个R代码文件life.R（注释文档后加一个NULL，再加一个空白行）
#' @title
#'
#' @description
#' @docType data
#' @keywords datasets
#' @name life
#' @usage life
#' @format
#' @source
NULL
#R要求每一个包都有严格的结构化文档
#LaTeX（标记语言和排版系统）编写文档，roxygen2包来简化文档创建过程
#创建包步骤：
#①安装工具：安装roxygen2包，Rtools.exe和MIKTeX工具（Windows）
#②配置文件夹：R(oneway.R,print.R,summary.R,plot.R,npar.R,life.R),data(life.rda)
#③生成文档：roxygenize("npar")
#namespace命名空间控制函数的可视性
#④建立包：system("R CMD build npar") 或system("Rcmd INSTALL --build npar")建立二进制包
#⑤检查包（可选）：system("R CMD check npar")发布到CRAN需无任何警告
#⑥创建pdf手册（可选）：system("R CMD Rd2pdf npar")
#⑦本地安装包（可选）：system("R CMD INSTALL npar")
#⑧上传到CRAN（可选）------------------第22章 创建动态报告--------------------------------------
#四种类型的模板：R Markdown/ ODT/ DOCX/ LaTeX(出版级别文档，语法较复杂，学习曲线陡峭)
#1.R markdown
# ```{r echo=F,results='hide'}
# ````
#R代码块参数：
# echo：T/F是否输出源码
# results：asis/hide是否输出结果
# warning：T/F是否输出警告
# message：T/F是否输出参考信息
# error：T/F是否输出错误信息
# fig.width
# fig.height
#用Rstudio创建和处理R markdown
#a. File——new File——R markdown 生成骨架文件进行编辑
#b. 再点击Knit——render(to Html/PDF/Word)
#PS: 生成PDF文档需要安装MiKTeX（windows）
#----------------第23章 使用lattice进行高级绘图---------------------------
#特长：擅长网格图形，展示变量的分布或变量间的关系，每幅图代表一个变量或各个水平
library(lattice)
histogram(~height | voice.part,data=singer)
#height是独立变量，voice.part是调节变量
#y~x|A*B
#x为主要变量，AB为调节变量。小写字母代表数值型变量，大写字母代表分类变量（因子）
#lattice图类型及相应的函数：
contourplot()
levelplot()
cloud()
wireframe()
barchart()
bwplot()
dotplot()
histogram()
densityplot()
parallelplot()
xyplot()
splom()
stripplot()
#例子：
library(lattice)
attach(mtcars)
gear <- factor(gear,levels = c(3,4,5),labels = c("3 gears","4 gears","5 gears"))
cyl <- factor(cyl,levels = c(4,6,8),labels = c("4 cylinders","6 cylinders","8 cylinders"))
head(mtcars)
densityplot(~mpg)
densityplot(~mpg|cyl)
bwplot(cyl~mpg|cyl)
xyplot(mpg~wt|cyl*gear)
cloud(mpg~wt*qsec|cyl)
dotplot(cyl~mpg|gear)
splom(mtcars[c(1,3,4,5,6)])
detach(mtcars)
#图形参数：
mygraph <- densityplot(~height|voice.part,data=singer);mygraph
update(mygraph,col="red",pch=16,cex=0.8,jitter=0/05,lwd=2) #update函数调整图形对象
#lattice绘图的一个强大特征是可以增加调节变量，通常调节变量是因子
#若是连续变量，可先转换为shingle结构，再作为一个调节变量
displacement <- equal.count(mtcars$disp,number=3,overlap=0) #将连续变量disp分成3个无重叠的间隔
xyplot(mpg~wt|displacement,data=mtcars,
layout=c(3,1),aspect = 1.5) #布局和纵横比（高/宽）
#面板函数
#每一个绘图函数都对应一个绘图面板panel.graph_function，如上图也可写成：
xyplot(mpg~wt|displacement,data = mtcars,panel=panel.xyplot)
#因此我们可将lattice默认的50多个默认函数组合成自定义函数，便类似于ggplot2的图层
mypanel <- function(x,y){
panel.xyplot(x,y,pch=19) #散点图
panel.rug(x,y) #地毯图
panel.grid(h=-1,v=-1) #网格线
panel.lmline(x,y,col="red",lwd=1,lty=2) #回归曲线
}
xyplot(mpg~wt|displacement,data=mtcars,
layout=c(3,1),aspect = 1.5,
panel = mypanel)
#例2
mtcars$transmission <- factor(mtcars$am,levels = c(0,1),labels = c("automatic","manual"))
panel.smoother <- function(x,y){
panel.grid(h=-1,v=-1)
panel.xyplot(x,y)
panel.loess(x,y) #非参拟合曲线
panel.abline(h=mean(y),lwd=2,lty=2,col="darkgrey") #水平参考线
}
xyplot(mpg~disp|transmission,data=mtcars,
scales = list(cex=0.8,col="red"),
panel=panel.smoother)
#分组变量
densityplot(~mpg,data=mtcars,auto.key = T, #图例放上方
group=transmission) #添加分组变量每个水平的图
#一个带有分组和调节变量以及自定义图例的例子：
head(CO2)
colors <- "darkgrey"
symbols <- c(1:12)
linetype <- c(1:3)
#自定义图例
key.species <- list(title="Plant",
space="right",
text=list(levels(CO2$Plant)),
points=list(symbols,col=colors))
#Plant是分组变量，Type和Treatment是调节变量
xyplot(uptake~conc|Type*Treatment,data=CO2,
group=Plant,tpye="o",
pch=symbols,col=colors,lty=linetype,
key=key.species)
show.settings(xyplot)
show.settings() #查看图形参数默认设置(图形化)
mysettings <- trellis.par.get() #得到图形默认设置
names(mysettings)#查看图形参数列表成分
mysettings$grid.pars
mysettings$superpose.symbol #参数中具体某一成分的默认值
mysettings$superpose.symbol$pch <- c(1:10) #改变默认值
trellis.par.set(mysettings) #更改参数设置
show.settings(mysettings) #查看改动
#自定义图形条带：strip.custom
#页面布局：split=c(1,2,1,2)-----------------------附录------------------------------
#R GUI（如R Commander）与SAS和SPSS相比，不丰富和不成熟。
#Shiny实现web互动
#自定义启动环境：通过修改Rprofile.site（站点初始化文件）和.Rprofile（目录初始化文件），R启动时会执行其中代码
#R启动步骤：加载R_HOME/etc/Rprofile.site文件——> 当前目录寻找.Rprofile，若没找到该文件，则到用户主目录HOME中去找
Sys.getenv("R_HOME") #R_HOME环境变量位置
Sys.getenv("HOME") #用户主目录HOME
getwd() #当前工作目录
#可以在这两个文件中放入两个特殊函数：.First函数（R会话开始时执行）和.Last函数（会话结束时执行）
#Rprofile.site文件自定义示例：
#常用选项设置
options(papersize = "a4")
options(editor = "notepad")
options(tab.width=2)
options(width = 130)
options(digits = 4)
options(stringsAsFactors = FALSE)
options(show.signif.stars = FALSE)
grDevices::windows.options(record=TRUE)
#设置R交互提示符
options(prompt = ">")
options(continue = "+")
#设置本地库路径
.libPaths("D:/my_R-library") #在R目录外创建一个本地库
#设置默认的CRAN镜像
local({r <- getOption("repos")
r["CRAN"] <- "http://cran.cae.edu/"
options(repos=r)
})
#启动函数
.First <- function(){
library(lattice)
library(ggplot2)
source("E:/mydir/myfunctions.R")
cat("\nWelcome at",date(),"\n")
}
#会话结束函数
.Last <- function(){ #可执行某些清理安装，如保存历史记录、保存程序输出、保存数据等
cat("\nGoodbye at",date(),"\n")
}
#更多方法：
help("Startup")
##2.更新R
#升级R自身较为复杂，主要是升级完后很多自定义选项和扩展包会丢失
#手动将旧版R中包批量导入新R中的方法：
oldip <- installed.packages()[,1];oldip #保存旧版已安装的包
save(oldip,file="E:/installedPackages.Rdata") #保存已安装的包到非R目录下
load("E:/installedPackages.Rdata") #打开新版R后，载入旧版R包
newip <- installed.packages()[,1]
for(i in setdiff(oldip,newip)){
install.packages(i)
}
#但此方法只能安装CRAN上的包，如Bioconductor上的R包则只能重新安装
source("http://biocondutor.org/biocLite.R")
biocLite("globaltest")
biocLite("Biobase")手机扫一扫
移动阅读更方便
你可能感兴趣的文章