R　COOKBOOK　学习笔记

 
 help.start()
 args("lm") #快速获得函数的参数
 example("lm")
 help.search("lm")
 help.search("plot",package="party")
 help(package="party")
 vignette(package="party")  # 查看附加文档列表
 RSiteSearch("tree") #R网站查询
 
 print() cat() #print 一次只能打印一个对象 而cat 不是， 但cat 不能显示复合的对象， 如矩阵等
 x<<-3  #强制赋值给全区变量
 ls.str() #在ls的基础上带了对象信息
 rm(x)
 rm(list=ls())
 seq(from=1,to=5,by=2)
 seq(from=0,to=100,length.out=5)
 fib[-(1:3)]
 v[v>median(v)]
 v[(v<quantile(v,0.05))|(v>quantile(v,0.95))]
 v[!is.na(v)&is.NULL(v)]
 years["Carter"] #可使用名称引用向量中的元素
 (w-mean(w))/sd(w) #z~score
 
 require("party") #便于脚本的编写， 在载入失败时返回F
 library("party") # 失败是错误信息  此两函数相似
 
 getwd()
 setwd()
 save.image()
 
 .last.value #特殊变量，储存最近计算出的表达式值
 search() #查看当前R中载入的包
 
 getAnywhere("lda") #查找函数在哪个包里面
 MASS:::lda(x) #使用此函数
 data()
 data(Cars93,package="MASS") #载入包中的数据集
 library() #不含参数时列出已经安装的R包列表
 installed.packages()[,c("Package","Version")] # 列出安装的包
 source("hello.R",echo=T) #显示R提示符的命令代码
 
 Sys.getenv("SHELL")
 Sys.getenv("DISPLAY")
 Sys.getenv("R_HOME") #找到R的主目录
 Rprofile #设定初始参数，包括初始载入包等
 
 print(pi,digits=4)
 cat(format(pi,digits=4),"\n")
 
 sink("script_output.txt")  #输出结果到文件
 source("script.R")
 sink() # 关闭sink
 
 cat(result,file="analysisRepart.out",append=T) 
 con<-file("analysisReport.out","w") #文件连接
 cat(data,file=con)
 cat(results,file=con)
 close(con)
 
 list.files()
 records<-read.fwf("fixed-width.txt",widths=c(10,10,4,-1,4)，col.names=c("Last","First","Born","Died") #读取固定宽度的文件
 dfrm<-read.table("statisticians.txt",stringAsFactor=F)
 tbl<-read.csv("table-data.csv",as.is=F)  #不转换为因子
 write(tbl,file="table-data.csv",row.names=T)
 
 library(xml)
url<-'http://www.example.com/data/table.html'
tbls<-readHTMLTable(url)
tbl<-readHTMLTable(url,which=3)


read.table, read.csv 面向行读取数据
readLines, scan 能面向行甚至每一个单元读取
lines<-readLines("input.txt",n=10)
singles<-scan("singles.txt",what=numeric(0))
triples<-scan("triples.txt",what=list(character(0),numeric(0),numeric(0))


world.series<-scan("http://lib.stat.cmu.edu/datasets/wseries", skip=35,nlines=23,what=list(year=integer(0),pattern=character(0)),)
library(RMySQL)
con<-dbConnect(MySQL(),user="userid",password="pswd",host="hostname",client.flag=CLIENT_MULTI_RESULT)
sel<-"select * from SurveyResults where city="Chicago"
rows<-dbGetQuery(con,sql)
if(dbMoreResults(con)) dbNextResults(con)
dbDisconnect(con)


save(MyData, file="MyData.Rdata")
load("MyData.RData)
dput, dump #采用ASCII保存


dim(A)<-c(2,3)
B<-list(1,2,3,4,5,6)
dim(B)<-c(2,3)


append(1:10,99,after=5) #向量合并比较慢
cbind(1:6,1:3) # 向量循环相加


unstack() stack() # 数据融合， 数据框按照向量合并成两列，然后又转化为数据框
years[c("Carter","Clinton")]<-NULL   #移除列表元素
unlist() #列表转换为向量
lst[sapply(lst,is.null)]<-NULL # 从列表中移除空值？？， 注意Sapply应用于列表的每一个元素， 注意和years[c("Carter","Clinton")]<-NULL 区分
 
 dim(v)<-c(2,3) #向量直接变为矩阵
 solve(A) # 矩阵求逆
 
 
 rownames(mat)<-c("job","department")
 colnames(mat)<-c("job","department")
 
 row<-mat[1,drop=FALSE] #注意drop
 as.data.frame(lst)
 do.call(rbind.obs)
 
  f<-list(list(c(1,2,3),c(2,3,4)),list(c(6,7,8)))
  Map(as.data.frame,f)  # Map 把”列表的列表“转为数据框
  do.call(rbind,Map(as.data.frame,obs)
  
  N<-100000
  dfrm<-data.frame(dosage=numeric(N),lab=character(N))  #预先分配数据类型和行数
  subset(dfrm,select=c(predictor,response),subset=c(response>0)) # 注意应用，选择了列，选择了行
  subset(Car93,select=c(Manfufacturer,Model),subset=c(MPG.highway>median(MPG.highway)))
  subset(patient.data,select=c(-patient.id,-dosage))
  rbind(vec),cbind(vec) #建立一行或一列矩阵
  
  
  
  edit(dfrm)
  fix(dfrm)
  
  rbind() cbind() #也可进行数据框的操作
  merge(born,died,by="name") #数据框的连接  参数选择，类似SQL join 
  with(suburbs,(pop-mean(pop))/sd(pop)) #更容易的使用数据框列变量
  
  attach(suburbs)
  suburbs$pop<0
  pop ##插入一个数据框到搜索列表会插入一个临时副本，意味着对原数据框的改变被隐藏了，容易引起混淆
  
  
  clean<-na.omit(dfrm)   #数据框移除na， 列表不可以
  
  
  lapply(), sapply() # lapply 的返回结果是列表，而sapply的返回结果是向量，或者矩阵
 cor<-sapply(pred,cor,y=resp) # 注意第三个参数是第二个参数的参数
 tapply(x,f,fun) #将函数应用于列分组 f 为因子
 split()
 by(rials,trials$sex,summary) #把函数应用于行分组
 mapply(f, wec1,wec2,..wecn) #把函数应用于向量的每一个数
 
 nchar("James") #获取字符串长度 length 返回的是向量的长度
 s<-c("James","is","Good","emploee")
 substr(s,2,3)   #抽取字符
 
  strsplit(path,"/") #第二个参数可以是正则表达式 结果是列表
  sub(old,new,string) , gsub(old,new,string) # 字符代替，gsub代替所有实例， sub代替第一个实例
  print(), cat() #print 能显示出转移字符等特殊字符， 而cat不能，所以用print 打印特殊字符。
  
  f<-c("T1","T2","T3")
  d<-c("U1","U2","U3")
  m<-outer(f,d,paste,sep="-")   # outer本用于计算矩阵外积， 但可以变换函数， 可以生成卡迪尔积
  m[!lower.tri(m)] #上面产生的矩阵用下三角移除重复
  Sys.Date()
  as.Date("12/31/2010",format="%m/%d/%Y") # 注意大小写代表不同含义
  format(sys.Date())
  as.character(sys.Date())
  ISOdate(2012,2,29)
  as.date(ISOdate(2012,2,29))
  d<-as.Date("2010-3-15")
  p<-as.POSIXlt(d)   # 分解日期 年 月 日 等。。
  seq(from=s,by="month",length.out=12) #创建时间序列
  
  choose(5,2) #计算组合数
  combn(1:5,3) #生成组合
  runif(10) #均匀分布
  sample(c(F,T),20,replace=T,prob=c(0.2,0.8))
  d,p,r #d 开头概率密度函数，也称分布函数， p 开头累积概率函数
  pbinom(7,size=10,prob=0.5,lower.tail =F) # 右尾概论
  diff(pbinom(c(3,7),size=10,prob=0.5))  # 计算向量差， 计算两个概率之差
  qnorm(c(0.025,0.975)) # 分位数
  
  x<-seq(from=0,to=6,length.out=100)
  ylim<-c(0,0.6)
  plot(x,dunif(x,min=2,max=4),main="Uniform",type='l',ylim=ylim)
  abline(h=0.45,v=2)
  ploygon(region.x, region.y,density=10) #增加阴影区域
  mean(x>1) # 用这样的方式计算频数， 先转化成1，0 组成的向量，然后计算均值，最后就是发生的频数
  table(job, age) ，xtabs #列联表
  quantile(vec,c(.05.95)) #分位数
  z<-scale(f) # Z分数
  x<-rnorm(50,mean=100,sd=15)
  t.test(x,mu=95)
  t.test(x,conf.level=0.99)
  wilcox.test(x,conf.int=T) # 中位数置信区间
  prop.test(x,n,p) #样本比例检验
  runs.test(as.factor(s))  #游离检验
  t.test(x,y) #比较两样本均值
  t.test(x,y,paired=T) 
  wilcox.test(x,y,paired=T)  #
  wilcox.test(fav,unfav,paired=T) #
  cor.test(x,y,method="Spearman")
  success<-c(14,10)
  trials<-c(38,40)
  prop.test(success,trials)
  prop.test(x,n,p)
  pariwise.t.test(x,f)
  ks.test(x,y)
  plot(x,main="Forecast Result",xlab="Month",ylab="Production",col=c("red","black","green"))
  
  plot(x,main="the title",xlab="X-axis Label",ylab="Y-axis Label")
  plot(x,ann=F) # ann 表示先不绘制注释内容，然后调用title添加
  title(main="The Title",xlab="X Axis Label",ylab="Y Axis Label")  #增加内容
  grid() #网格
  points(13,21)
  with(iris,plot(Petal.Length,Petal.Width,pch=as.integer(Species)))
  f<-factor(iris$Species)
  legend(1.5,2.4,as.character(levels(f)),pch=1:length(levels(f)))
  abline(v=6)
  t<-lm(x[,1]~x[,2])
  abline(t)  #增加回归线
  plot(iris[,1:4]) # 多个变量，两两作图
   data(Cars93,package="MASS")
   coplot(Horsepower~MPG.city|Origin,data=Cars93) # 条件化图 condition plot |后面是因子
   heights<-tapply(airquality$Temp,airquality$Month,mean)
   barplot(heights,main="Mean Temp.by Month",
   names.arg=c("May","Jun","Jul","Aug","Sep"),ylab="Temp (deg.F)")  #条形图
   barplot2(heights,plot.ci=T,ci.l=lower,ci.u=upper) #增加置信区间的条形图， 置信区间需要单独算出然后带回
   rel.hts<-(heights-min(heights))/(max(heights)-min(heights))
   grays<-gray(1-rel.hts) #灰度  rainbow
   barplot(heights,col=grays,ylim=c(50,90),xpd=F,main="Mean Temp.By Month",names.arg=c("May","Jun","Jul","Aug","Sep"),ylab="Temp (deg.F)")
   plot(x,y,type="l",lwd=2)
   lines(x,y.indeps,col="yellow")
   
  xlim<-range(c(x1,x2))
  ylim<-range(c(y1,y2))
  plot(x1,y1,type='l'xlim=xlim,ylim=ylim)
  lines(x2,y2,lty="dashed")
  data(UScereal,package="MASS")
  boxplot(sugars~shelf,data=UScereal)
  lines(density(samp))
  plot(table(x),type='h',lwd=5,ylab="Freq")
  qqnorm(x)
  qqline(x)
 qqnorm(log(Cars93$Price),main="QQ plot: price")
 qqline(log(Cars93$Price))
 
 RATE<-1/10
 plot(qexp(ppoints(y),rate=RATE),sort(y))
 abline(a=0,b=1)
 colors<-ifelse (x>=5,"black","gray")
 plot(x, type='h',lwd=3,col=colors)
 f<-function(x) exp(-abs(x))*sin(2*pi*x)
 curve(f,-5,+5,main="Dampened Sine Wave")
 
 par(mfrow=c(2,3)) # 图像划分为2*3的矩阵
 par(mfcol=c(2,3)) # 按列增加图形
 win.graph(width=7.0,height=5.5) #打开另外一个图形窗口
  dev.set(2) #  设置活动窗口
 savePlot(filename="filename.ext",type=type)
 png("Mplot.png",width=648,height=432)
  curve(sin,-3,+3)
  dev.off()
  par(lwd=2)  # 许多设置参数， cex,col bg.
  
  lm(y~x)
  lm(y~u+v+w)
  lm(y~u+v+w,data=dfrm)
  
  m<-lm(y~x)
  coefficients(m)
  confint(m)
  deviance(m)
  fitted(m)  ##拟合y值的向量
  vcov(m) #协方差矩阵
   anova(m)
  lm(y~x+0) ##截距为零的线性方程
  lm(y~u*v) # 交互项的线性方程 包括所有可能的交 相当于y~u+v+u*v
  lm(y~u+v+w+u:v:w) # 只包含了uvw的交互，相当于y~u+v+w+uvw 没有uv,uw,wv 
  lm(y~(u+v+w)^2) # 二阶互交， 也可三阶， *，： 符合分配率， x*(u+v+z),x:(u+v+z)
  
  full.model<-lm(y~x1+x2+x3+x4)
  reduced.model<-step(full.model,direction="backward")
  
  min.model<-lm(y~1)
  fwd.model<-step(min.model,direction="forward",scope=(~x1+x2+x3+x4))
  fwd.model<-step(min.model,direction="forward",scope=(~x1+x2+x3+x4),trace=0) #只显示选择后的
  
  lm(y~x,subset=1:500)
  lm(y~x,subset=1:floor(length(x)/2))   # 选择子集进行回归
  lm(y~x,subset=(lab=="NJ"))
  m=lm(y~u+I(u^2))  #I(...) 计算表达式然后进行回归
  predict(m,newdata=data.frame(u=13.4))
  lm(y~poly(x,3,raw=TRUE)) # 多项式回归， raw为T
  lm(y~x+I(x^2)+I(x^3))
  lm(log(y)~x) 
  lm(sqrt(y)~month)
  lm(y~sqrt(x))
  lm(log(y)~log(x))
  library(MASS)
  m<-lm(y~x)
  boxcox(m) # 确定指数，然后进行回归。当线性无法满足了
  which.max(bc$y)
  lambda<-bc$x[which.max(bc$y)]
  z<-y^lambda
  ms<-lm(z~x)
  library(car)
  outlier.test(m) # 利用car包中的outlier来检测离群值
  influence.measures(m) # 识别有影响的观察值，异常？
  library(lmtest)
   dwtest(m)  # Durbin-Watson检查残差自相关
   dwtest(m,alternative="two.side")
   acf(m)  #检查残差自相关
   
   m<-lm(y~u+v+w)
   preds<-data.frame(u=3.1,v=4.0,w=5.5)
   predict(m,newdata=preds) 
   predict(m,newdata=preds,interval="prediction") 
   
   oneway.test(x~f)
   m<-aov(x~f)
   summary(m)
   
  library(faraway)
  data(rats)
  interaction.plot(rats$poison,rats$treat,rats$time)  #互交图
  
  m<-aov(x~f)  ##只进行方差分析验证均值不一样
  TukeyHSD(m)    ### 显示出那俩组数值的均值不一样。
  plot(TukeyHSD)
  
  kruskal.test(x~f) #不满足正态分布假设的费参数检验，稳健的检验中位数是否相同
  anova(m1,m2 ) ##比较两个模型是否有显著的区别
  
  tail(x) #查看最后几行
  option(width=120) #设置输出窗口的宽度
  
  x<-rnorm(1000)
  breaks<-c(-3,-2,-1,0,1,2,3)
  labels<-c("Bottom","Low","Neg","Pos","Hight","Top")
  f<-cut(x,breaks,labels) # 分箱
  summary(f)
  
  vec<-c(1,2,3,4,5)
  match(3,vec)
  
  v[seq_along(v)%%n==0]  ##每隔n个值取一个变量
  
  pmin()
  pmax()
  
  sides<-factor(c("Heads","Tails"))
  faces<-factor(c("1 pip",paste(2:6,"pips")))
  expand.grid(faces,sides)  ##生成多个因子的组合
  
  as.vector(as.matrix(dfrm))
  mean(as.matrix(daily.prod))
  
  dfrm<-dfrm[order(dfrm$key1,dfrm$key2),] # 对数据框两列排序
  
  m<-lm(resp~pred)
  slop<-coef(m)[2]
  str(slop)
  attributes(slope)<-NULL  # 移除变量属性
  attr(slop,"name")<NULL   # 移除变量属性
  
  class() #看是什么对象
  mode() #移除对象特性看内不， 它把class的特性去掉了
  str()  #看内部结构
  
  system.time(sum(rnorm(1000000)))
  warnings()
  suppressWarnings()
  
  numbers<-list(1,3,5,7,9)
  mean(unlist(numbers))
  
  lists<-list(col1=list(7,8,9),col2=list(70,80,90),col3=list(700,800,900))
  cbind(lists) # 不能很好的转化为矩阵
  do.call(cbind,lists)  ##转化为了矩阵
  cbind(lists[[1]],lists[[2]],lists[[3]])   ##转化为了矩阵 
  
  '%+-%'<-function(x,margin) x+c(-1,+1)*margin  # 定义二元运算符
  
  f<-function(x)3*x^4-2*x^3+3*x^2-4*x+5
  optimize(f,lower=-20,upper=20)  # 求单参数函数最大值和最小值
  optimize(f,lower=-20,upper=20,maximum=TRUE)
  
   
  f<-function(v){ a<-v[1];b<-v[2]
                  sum(abs(z-((x+a)^b)))
 }
    optim(c(1,1),f) # 求多参数函数最大值和最小值 c(1,1)为开始搜索点
optim(c(1,1),f control=list(fnscale=-1))   
eigen() #计算特征值和特征向量

  r<-prcomp(~x+y+z)  # 主成分分析
       summary(r)
  print(r)

   r<-prcomp(~x+y)
       slope<-r$rotation[2,1]/r$rotation[1,1]  #简单正交回归 主成分载荷计算斜率
  intercept<-r$center[2]-slope*r$center[1] # 用斜率在计算tie
  
  hc<-hclust(d)
  clust<-cutree(hc,k=3)
       head(clust,20)
       plot(x~factor(means),main="Original Clusters",xlab="Cluster Mean")
  
      data(pima,package="faraway")
      b<-factor(pima$test)
 m<-glm(b~diastolic+bmi,family=binomial, data=pima)  # 逻辑回归
 
 m.red<-glm(b~bmi,family=binomial,data=pima)
      newdata<-data.frame(bim=32.0)    #以data.frame存放数据
      predict(m.red,type="response",newdata=newdata)
 
 
 
 stat<-function(data,indices){
      r<-prcomp(~x+y,data=data,subset=indices)
      slope<-r$rotation[2,1]/r$rotation[1,1]
      return(slope)
       }
     boot.data<-data.frame(x=x,y=y)
     library(boot)
     reps<-boot(boot.data,stat,R=999) # 自助法抽样
     boot.ci(reps,type=c("perc","bca")) # 自助法抽样并求置信区间
factanal(diffs,factors=2) #因子分析
 
library(zoo) #时间序列包
library(xts)  #时间序列包
    prices<-c(132.45,130.85,130.00,129.55,130.85)
    dates<-as.Date(c("2010-01-04","2010-01-05","2010-01-06","2010-01-07","2010-01-08"))
    ibm.daily<-zoo(prices,dates)  #创建时间序列对象，prices可以是数据框， 是多元时间向量
    print(ibm.daily)
plot(ibm.daily2,screens=c(1,2))
ibm.daily[as.Date('2010-01-05')]
dates<-seq(as.Date('2010-01-04'),as.Date('2010-01-08'),by=2)
ibm.daily[dates]
filed.cpi<-na.locf(merge(cpi,empty,all=T)) #合并时间序列并补全缺失值
lag(ibm.daily,k=+1,na.pad=T) #时间序列的滞后
diff(ibm.daily) #计算差分
diff(ibm.daily,lag=12) 
log(ibm.daily)
ma<-rollmean(ts,k, align="right") #计算移动平均
apply.daily(ts,f) #日历范围内使用函数





vignette("zoo") # 查看包文档
plot(sqrt(251)*apply.monthly(as.xts(diff(log(ibm))),sd),main="IBM:Monthly Volatility",ylab="Std dev,annualized")
sds<-rollapply(ts,30,sd,by=30,align="right") #滚动函数
acf(ts) #自相关
box.test(ts) #检验时间序列自相关性
pacf(ts) #偏自相关
ccf(ts1,ts2) #交叉相关函数， 滞后的相关性
m<-lm(coredata(ts)~index(ts))
detr<-zoo(resid(m),index(ts)) #利用线性回归找到趋势， 然后剔除趋势
library(forecast)
auto.arima(x)
arima(x,order=c(p,d,q)
m<-ariam(x,order=c(2,1,2),fixed=c(0,NA,0,NA) #剔除不显著的参数
tsdig(m) #对模型进行检验
predict(m,n.ahead=10) #预测
library(tseries)
adf.test(coredata(ts))  #均值回归检验
library(fUnitRoots)
adfTest(coredata(ts1),type="nc")
library(KernSmooth)
t<-seq(from=-10,to=10,length.out=201)
    noise<-rnorm(201)
   y<-sin(t)+noise
    gridsize<-length(y)
   bw<-dpill(t,y,gridsize=gridsize)
   lp<-locpoly(x=t,y=y,bandwidth=bw,gridsize=gridsize) ##对数据平滑处理
   smooth<-lp$y