看了r导论和r_for_beginner统计建模与r软件，笔记大致按照简介->运算符->数据结构->绘图->假设检验->回归分析顺序
r_intro.pdf
chapter1 绪论 #讲R的运行方式（对话式命令行）和用途
R 是开发新的交互式数据分析方法一个非常好的工具。
部分的统计功能是整合在 R 环境的底层，但是大多数功能则以包 的形式提供。大约有25个包和 R 同时发布（被称为\标准" 和\推荐" 包），更多的包可以通过网上或其
他地方的CRAN 社区(http://CRAN.R-project.org) 得到。
help(function) or ?function
基本命令要么是表达式（expressions）要么就是赋值（assignments）。
命令可以被(;)隔开，或者另起一行。基本命令可以通过大括弧(f和g) 放在一起构成一个复合表达式（compoundexpression）。
如果一条命令在一行结束的时候在语法上还不完整， R 会给出一个不同的提示符，默认是+
在每一次 R 会话结束的时候，你可以保存当前所有可用的对象。如果你想这样
做，这些对象将会写入当前目录下一个叫.RData10 的文件中，并且所有在这次会话中
用过的命令行都会被保存在一个叫.Rhistory 的文件中。
如果采用 R 做分析，你最好用相对独立的工作目录。
chapter2 简单的算术操作和向量运算
2.1 向量和赋值
其中，最简单的结构就是由一串有序数值构成的数值向量（vector）
>assign("x",c(10.4,5.6,3.1,6.4,21.7)) >c(10.4,5.6,3.1,6.4,21.7)->x
2.2 向量运算
+ - * / ^ sqrt max min length sum mean var ##square root
2.3 生产序列
c() 1:30 seq(form=2,to=30,length=15) seq(from=2,to=30,by=2) rep(x,times=5) rep(x,each=5)
2.4 逻辑向量
逻辑向量可以由条件式（conditions）产生 例如 >temp<-x>13
逻辑运算符 < <= > >= == != & |
2.5 缺损值
NA is.na(x) NaN in.nan(x)
2.6 字符向量
"" \t \n \b \\
函数paste() 可以有任意多的参数，并且把它们一个接一个连成字符串。
>labs<-paste(c("X","Y"),1:10,sep="")
>labs
c("X1","Y2","X3","Y4","X5","Y6","X7","Y8","X9","Y10")
2.7 索引向量
>y<-x[!is.na(x)]
>(x+1)[(!is.na(x))&x>0]->z
>x[1:10]
>y<-x[-(1:5)]
>fruit<-c(5,10,1,20)
>names(fruit)<-c("orange","banana","apple","peach")
>lunch<-fruit[c("apple","orange")]
其他数据结构：matrix, array, list, factor, data frame, function
第三章对象及它们的模式和属性
实数或复数向量，逻辑向量和字符串向量之类的对象属于\原子"（atomic）型的对象，因为它们的元素都是一样的类型或模式。
列表被认为是一种\递归"结构而不是原子结构,因为它们的元素可以以它们各自的方式单独列出。
3.1 内在属性：模式和长度
mode() numeric, logical, character, complex, raw
length()
3.2 改变对象长度
>length(alpha)<-3
3.4 对象的类
class() numeric logical, character, list, matrix, array, factor, data,frame
引入对象的类属性有利于面向对象风格的7 R 编程。比如说，如果一个对象属于"data.frame" 类，那么它将会以一种特定的方式显示，函数plot() 也会以特定的方式显示它的图形。
第四章有序因子和无序因子
因子（factor）是一个对等长的其他向量元素进行分类（分组）的向量对象。
4.1一个例子
>statef<-factor(state)
>levels(statef)
定义函数
>stderr<-function(x)sqrt(var(x)/length(x)) # sqrt=square root
作为一个练习，你可以计算一下州平均收入的95%信度区间。提示一下，你可能要再次使用tapply() 和能得到样本量的函数length()，以及能得到t-分布分位数的函数qt()(你需要参考一下 R 为t-检验设计的函数。)。
5.1数组
数组可以看作是带有多个下标类型相同的元素集合，如数值型。
维度向量（dimensionvector）是一个正整数向量。如果它的长度为k，那么该数组就是k-维的，例如矩阵是2-维数组。
定义数组 dim(z) matrix() array()
matrix(data = NA, nrow = 1, ncol = 1, byrow = FALSE,
dimnames = NULL)
array(data = NA, dim = length(data), dimnames = NULL)
5.2 数组索引
假定数组a的维数向量是c(3,4,2)，a[2,,] 是一个4 * 2 的数组，它的维度向量为c(4,2)。
5.3 索引数组
和用于下标位置的索引向量一样，可以根据索引数组去给数组中不规则的元素集合赋值或者将数组中特定的元素返回到一个向量中。举例：p36
5.4 array()函数
>Z<-array(data vector, dim vector)
假定向量h 有24个或更少的数值，那么命令
>Z<-array(h,dim=c(3,4,2))
矩阵
colnames() rownames() ncol() nrol()
t()
第六章列表和数据框
6.1 列表
R 的列表（list）是一个以对象的有序集合构成的对象。列表中包含的对象又称为它的分量（components）。
>Lst<-list(name="Fred",wife="Mary",no.children=3,child.ages=c(4,7,9))
这些分量则可以用Lst[[1]], Lst[[2]], Lst[[3]] 和Lst[[4]]独立访问。如果Lst[[4]] 是一个有下标的数组，那么Lst[[4]][1] 就是该数组的第一个元素。
Lst$name Lst[[1]] Lst[["name"]] 三者等价
分量名字可以简写2，原则是只要能很好的区分所有分量就行。因此Lst$coefficients可以简写为Lst$coe 而Lst$covariance 可以简写成Lst$cov。
6.2 构建和修改列表
>Lst<-list(name 1 =object 1 , ... , name m=object m)
列表扩充
>Lst[5]<-list(matrix=Mat)
列表连接
>list.ABC<-c(list.A,list.B,list.C)
6.3 数据框
数据框（dataframe）是一个属于"data.frame" 类的列表。
创建数据框
>accountants<-data.frame(home=statef,loot=incomes,shot=incomef)
>kalythos <- data.frame(x = c(20,35,45,55,70), n = rep(50,5),y = c(6,17,26,37,44))
attach() detach()
r_for_beginner
chapter3 数据操作
>read.table(file, header = FALSE, sep = "", quote = "\"’", dec = ".",
row.names, col.names, as.is = FALSE, na.strings = "NA",
colClasses = NA, nrows = -1,
skip = 0, check.names = TRUE, fill = !blank.lines.skip,
strip.white = FALSE, blank.lines.skip = TRUE,
comment.char = "#")
>read.csv(file, header = TRUE, sep = ",", quote="\"", dec=".", fill = TRUE, ...)
>scan(file = "", what = double(0), nmax = -1, n = -1, sep = "",
quote = if (sep=="\n") "" else "’\"", dec = ".",
skip = 0, nlines = 0, na.strings = "NA",
flush = FALSE, fill = FALSE, strip.white = FALSE, quiet = FALSE,
blank.lines.skip = TRUE, multi.line = TRUE, comment.char = "")
write.table(x, file = "", append = FALSE, quote = TRUE, sep = " ",
eol = "\n", na = "NA", dec = ".", row.names = TRUE,
col.names = TRUE, qmethod = c("escape", "double"))
>save(x, y, z,file= "xyz.RData")
函数save.image()是save(list =ls(all=TRUE),file=".RData")的一个简捷方式。
chapter4 R绘图
R提供非常多样的绘图功能。如想了解，可以输入:demo(graphics) 或者demo(persp)。
有两种绘图函数：高级绘图函数（high-level plotting functions）创建一个新的图形，低级绘图函数（low-level plotting functions）在现存的图形上添
加元素。
绘图参数（graphical parameters）控制绘图选项，可以使用缺省值或者用函数par修改。
4.1 管理绘图
在Unix/Linux 下，绘图窗口称为x11，而在Windows下称为windows。
4.2 高级绘图函数
plot()
pie()
boxplot()
hist()
bar()
qqplot()
共同参数
main= xlab= ylab= xlim= ylim= sub= type="p","l",etc. add=TURE axes=TURE
4.3 低级绘图函数
统计建模与R软件
chapter4 参数估计
4.1点估计
距法 和 最小二乘法
4.3区间估计
4.3.1一个正态总体的情况
均值u
interval_estimate1<-function(x,sigma=-1,alpha=0.05){
n<-length(x);xb<-mean(x)
if(sigma>=0){
tmp<-sigma/sqrt(n)*qnorm(1-alpha/2);df<-n
}
else{
tmp<-sd(x)/sqrt(n)*qt(1-alpha/2,n-1);df<-n-1
}
data.frame(mean=xb,df=df,a=xb-tmp,b=xb+tmp)
}
方差
interval_var1<-function(x,mu=Inf,alpha=0.05){
n<-length(x)
if(mu<Inf){
S2<-sum((x-mu)^2)/n;df<-n
}
else{
S2<-var(x);df<-n-1
}
a<-df*S2/qchisq(1-alpha/2,df)
b<-df*S2/qchisq(alpha/2,df)
data.frame(var=S2,df=df,a=a,b=b)
}
4.3.2两个正态总体的情况
均值差 t.test()
t.test(x,y=NULL,
alternative=c("two.sided","less","greater"),
mu=0,paired=FALSE,var.equal=FALSE,
conf.level=0.95,...)
配对数据
t.test(X-Y)
方差比
interval_var2<-function(x,y,
mu=c(Inf,Inf),alpha=0.05){
n1<-length(x);n2<-length(y)
if(all(mu<Inf)){
Sx2<-1/n1*sum((x-mu[1])^2);Sy2<-1/n2*sum((y-mu[2])^2)
df1<-n1;df2<-n2
}
else{
Sx2<-var(x);Sy2<-var(y);df1<-n1-1;df2<-n2-1
}
r<-Sx2/Sy2
a<-r/qf(1-alpha/2,df1,df2)
b<-r/qf(alpha/2,df1,df2)
data.frame(rate=r,df1=df1,df2=df2,a=a,b=b)
}
var.test()
4.3.3非正态总体的区间估计
u的估计
interval_estimate3<-function(x,sigma=-1,alpha=0.05){
n<-length(x);xb<-mean(x)
if(sigma>=0)
tmp<-sigma/sqrt(n)*qnorm(1-alpha/2)
else
tmp<-sd(x)/sqrt(n)*qnorm(1-alpha/2)
data.frame(mean=xb,a=xb-tmp,b=xb+tmp)
}
4.3.4单侧置信区间估计
chapter5 假设检验
5.1假设检验的基本概念
5.1.1基本概念
一部分称为拒绝域W，当样本落入拒绝域时，便拒绝H0；另一部分称为接受域，当样品本落入它时不拒绝H0。
5.1.2基本思想与步骤
构造拒绝域W的常用方法是寻找一个统计量g， g的大小可以反映对原假设H0有利或不利。
5.1.3假设检验的两类错误
a=P{否定H0|H0是真实的} b=P{接受H0|H0是错误的} pi=1-b=P{否定H0|H0是错误的}
5.2.1正态总体均值的假设检验
1单个总体的情况
双侧检验H0 : µ = µ0,H1 : µ = µ0
|Z|≥ Zα/2,
|T|≥ t 2(n − 1),
单侧检验H0 : µ ≤ µ0,H1 : µ>µ0 或( H0 : µ ≥ µ0,H1 : µ<µ0),
T ≥ tα(n − 1)(( T ≤−tα(n − 1))
Z ≥ Zα 或 ( Z ≤−Zα).
求p值
P_value<-function(cdf,x,paramet=numeric(0),side=0){
n<-length(paramet)
P<-switch(n+1,
cdf(x),
cdf(x,paramet),
cdf(x,paramet[1],paramet[2]),
cdf(x,paramet[1],paramet[2],paramet[3])
)
if(side<0)P
elseif(side>0)1-P
else
if(P<1/2)2*P
else2*(1-P)
}
t检验程序
mean.test1<-function(x,mu=0,sigma=-1,side=0){
source("P_value.R")
n<-length(x);xb<-mean(x)
if(sigma>0){
z<-(xb-mu)/(sigma/sqrt(n))
P<-P_value(pnorm,z,side=side)
data.frame(mean=xb,df=n,Z=z,P_value=P)
}
else{
t<-(xb-mu)/(sd(x)/sqrt(n))
P<-P_value(pt,t,paramet=n-1,side=side)
data.frame(mean=xb,df=n-1,T=t,P_value=P)
}
}
2.两个总体的情况
方差已知
方差未知但相等
方差未知且不等
mean.test2<-function(x,y,
sigma=c(-1,-1),var.equal=FALSE,side=0){
source("P_value.R")
n1<-length(x);n2<-length(y)
xb<-mean(x);yb<-mean(y)
if(all(sigma>0)){
z<-(xb-yb)/sqrt(sigma[1]^2/n1+sigma[2]^2/n2)
P<-P_value(pnorm,z,side=side)
data.frame(mean=xb-yb,df=n1+n2,Z=z,P_value=P)
}
else{
if(var.equal==TRUE){
Sw<-sqrt(((n1-1)*var(x)+(n2-1)*var(y))/(n1+n2-2))
t<-(xb-yb)/(Sw*sqrt(1/n1+1/n2))
nu<-n1+n2-2
}
else{
S1<-var(x);S2<-var(y)
nu<-(S1/n1+S2/n2)^2/(S1^2/n1^2/(n1-1)+S2^2/n2^2/(n2-1))
t<-(xb-yb)/sqrt(S1/n1+S2/n2)
}
P<-P_value(pt,t,paramet=nu,side=side)
data.frame(mean=xb-yb,df=nu,T=t,P_value=P)
}
}
成对数据
>t.test(X-Y,alternative="less")
5.2.1正态总体方差的假设检验
1 一个总体的情况
var.test1<-function(x,sigma2=1,mu=Inf,side=0){
source("P_value.R")
n<-length(x)
if(mu<Inf){
S2<-sum((x-mu)^2)/n;df=n
}
else{
S2<-var(x);df=n-1
}
chi2<-df*S2/sigma2;
P<-P_value(pchisq,chi2,paramet=df,side=side)
data.frame(var=S2,df=df,chisq2=chi2,P_value=P)
}
2 二个总体的情况
var.test2<-function(x,y,mu=c(Inf,Inf),side=0){
source("P_value.R")
n1<-length(x);n2<-length(y)
if(all(mu<Inf)){
Sx2<-sum((x-mu[1])^2)/n1;Sy2<-sum((y-mu[2])^2)/n2
df1=n1;df2=n2
}
else{
Sx2<-var(x);Sy2<-var(y);df1=n1-1;df2=n2-1
}
r<-Sx2/Sy2
P<-P_value(pf,r,paramet=c(df1,df2),side=side)
data.frame(rate=r,df1=df1,df2=df2,F=r,P_value=P)
}
var.test(x,y,ratio=1,
alternative=c("two.sided","less","greater"),
conf.level=0.95,...)
二项分布总体的假设检验
binom.test(x, n, p = 0.5,
alternative = c("two.sided", "less", "greater"),
conf.level = 0.95)
5.2若干重要的非参数检验
这种不家丁总体分布的具体形式，尽量从数据（或样本）本身来获得所需要的信息的统计方法称为被阐述方法。
5.3.1pearson拟合优度卡方检验
K=sum(ni-npi)^2/npi
K ~ X^2(m-1)
可将P值称为所得数据域元假设的拟合优度，P值越大，支持元假设的证据就越强，给定一个显著性水平a， 当P值<a，就拒绝原假设。
chisq.test(x,y=NULL,correct=TRUE,
p=rep(1/length(x),length(x)),rescale.p=FALSE,
simulate.p.value=FALSE,B=2000)
检验正态分布
X<-scan()
25 45 50 54 55 61 64 68 72 75 75
78 79 81 83 84 84 84 85 86 86 86
87 89 89 89 90 91 91 92 100
A<-table(cut(X,br=c(0,69,79,89,100)))
p<-pnorm(c(70,80,90,100),mean(X),sd(X))
p<-c(p[1],p[2]-p[1],p[3]-p[2],1-p[3])
chisq.test(A,p=p)
检验poison分布
X<-0:6;Y<-c(7,10,12,8,3,2,0)
q<-ppois(X,mean(rep(X,Y)));n<-length(Y)
p[1]<-q[1];p[n]<-1-q[n-1]
for(iin2:(n-1))
p<-q-q[i-1]
chisq.test(Y,p=p)
k-s检验
ks.test(X,"pexp",1/1500)
ks.test(X,Y)
6 回归分析
6.1一元线性回归
6.1.1数学模型
yi = β0 + β1xi + εi i =1, 2,...n,
E(εi)=0, var(εi)= σ2 ,i =1, 2,...n.
6.1.2
回归参数的估计
最小二乘法
6.1.3
显著性检验
t法 F检验 相关系数
R中相关函数
lm() summary() anova() predict()
lm(formula, data, subset, weights, na.action,
method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
singular.ok = TRUE, contrasts = NULL, offset, ...) ##data可为data.frame或list,包含formula的变量。
>lm.sol<-lm(y~1+x)
>summary(lm.sol)
Call:
lm(formula=y~1+x)
Residuals:
Min1QMedian3QMax
-2.0431-0.70560.16940.66332.2653
Coefficients:
EstimateStd.ErrortvaluePr(>|t|)
(Intercept)28.4931.58018.045.88e-09***
x130.8359.68313.519.50e-08***
---
Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1
Residualstandarderror:1.319on10degreesoffreedom
MultipleR-Squared:0.9481,AdjustedR-squared:0.9429
F-statistic:182.6on1and10DF,p-value:9.505e-08
6.1.5预测
>new<-data.frame(x=0.16)
>lm.pred<-predict(lm.sol,new,interval="prediction",level=0.95)
>lm.pred
fitlwrupr
[1,]49.4263946.3662152.48657
6.1.7一个实例
forbes数据
lm.sol<-lm(log100~F,data=forbes)
>abline(lm.sol)
>y.res<-residuals(lm.sol);plot(y.res)
>text(12,y.res[12],labels=12,adj=1.2)
6.2R中与线性模型有关的函数
lm(formula,data,subset,weights,na.action,
method="qr",model=TRUE,x=FALSE,
y=FALSE,qr=TRUE,singular.ok=TRUE,
contrasts=NULL,offset,...)
(1) anova()
(2) coefficients()
(3) deviance()
(5) plot()
(6) predict()
(8) residuals()
(9)step()
(10) summary()
6.3多元线性回归分析
模型
Y = β0 + β1X1 + ··· + βpXp + ε
Y = Xβ + ε,
预测
>new<-data.frame(x1=80,x2=40)
>lm.pred<-predict(lm.sol,new,interval="prediction",level=0.95)
>lm.pred
fitlwrupr
[1,]123.9699117.2889130.6509
修正拟合模型
new.model<-update(old.model,new.formula)
fm5<-lm(y~x1+x2+x3+x4+x5,data=production)
fm6<-update(fm5,.~.+x6)
smf6<-update(fm6,sqrt(.)~.)
牙膏实例
>attach(toothpaste)
>plot(Y~X1);abline(lm(Y~X1))
>lm2.sol<-lm(Y~X2+I(X2^2))
>x<-seq(min(X2),max(X2),len=200)
>y<-predict(lm2.sol,data.frame(X2=x))
>plot(Y~X2);lines(x,y)
>lm.new<-update(lm.sol,.~.+I(X2^2))
>summary(lm.new)
Call:
lm(formula=Y~X1+X2+I(X2^2),data=toothpaste)
>lm2.new<-update(lm.new,.~.-X2)
>summary(lm2.new)
>lm3.new<-update(lm.new,.~.+X1*X2)
>summary(lm3.new)
6.4逐步回归
step()
step(object,scope,scale=0,
direction=c("both","backward","forward"),
trace=1,keep=NULL,steps=1000,k=2,...)
水泥实例
lm.sol<-lm(Y~X1+X2+X3+X4,data=cement)
>lm.step<-step(lm.sol)
>summary(lm.step)
add1(object,scope,...)
drop1(object,scope,...)
>drop1(lm.step)
>lm.opt<-lm(Y~X1+X2,data=cement);summary(lm.opt)
6.5回归诊断
6.5.1例图的重要性
6.5.2残差
普通残差 标准化残差 学生化残差
6.5.3残差图
1回归值与残差的残差图
>y.res<-resid(lm.sol);y.fit<-predict(lm.sol)
>plot(y.res~y.fit)
2残差的QQ图
>plot(model,2)
3变量与残差的残差图
>y.res<-resid(lm.sol)
>plot(y.res~x1);plot(y.res~x2)
plot(x,which=1:4,
caption=c("ResidualsvsFitted","NormalQ-Qplot",
"Scale-Locationplot","Cook’sdistanceplot"),
panel=points,
sub.caption=deparse(x$call),main="",
ask=prod(par("mfcol"))<length(which)&&dev.interactive(),
...,
id.n=3,labels.id=names(residuals(x)),cex.id=0.75)
6.5.4影响分析
帽子矩阵 hatvalues()
dffits()
cook统计量cooks.distance()
COVRATIO值 covratio()
influence.meature()
