---

R 语言实战（第二版）

part 5-2 技能拓展

---

----------第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")