基础统计与R作图.doc

一、数值变量的统计描述在R中没有一个直接的统计描述命令，可以把常见的指标都一起算出来。如概述中所说，R往往是提供一些比较基础的统计命令，需要使用者自己编程组合。1、 算术均数：mean(x)例如：c(12,13,14)-xmean(x) 将输出x的均数13。2、 中位数：median(x)3、 几何均数：我没有找到计算几何均数的命令，但可以通过对数函数来算几何均数。例如：log(x,10)-y 将x取自然对数mean(y)-z 计算y的算术平均数10z 将z取反对数，即得x的几何平均数4、 极差：max(x)-min(x)5、 样本方差：var(x)6、 样本标准差：sd(x)请注意都是“样本”。7、 变异系数：(sd(x)/mean(x)*100通过以上命令和程序，可以求出想要的描述指标。二、分类变量的统计描述1、 table(x) 该命令可以输出变量中各个值的绝对频数。2、 table(x/sum(x) 可输出相对数的频数。这个方法一般，还要把频数*百分比，不过我没有找到更好的方法。3、 table(x,y) 可输出绝对频数的交叉表 分类资料的描述，我还没有找到较好的命令。如有人知道，请指教。三、正态性检验可采用Shapiro-Wilk检验。程序如下：library(ctest)先加载ctest模块，许多检验方法都在此模块中。shapiro.test(x)可对x进行正态性检验。四、方差齐性检验程序如下：var.test(x,y)可对x和y进行方差齐性检验。五、统计图制作R是统计制图的高手。1、 直方图：hist例如：rnorm(1000)-gg产生样本量为1000的正态样本hist(gg)制图 结果如下：2、 饼图：pie例如：c(1,2,2,5)-x赋值y-x/sum(x)计算x的构成比names(y)xc(6,5,4,8,7,4,6,5,8,4,7,1,4,5)-yplot(x,y) 结果如下：4、 箱式图：boxplot例如：rnorm(10)-xboxplot(x)结果如下： R的制图功能非常强大，是R的亮点之一。一、变量数据直接输入1、数值变量如果有一个变量，变量名为x，其数据如下：12，15，46，23，15。命令语句如下：xx其中：x为变量名；为赋值符；c( )为向量建立函数，表示把括号中的数据建立为一个向量。以上命令语句建立了一个数值变量x。2、字符变量字符变量的建立与数值变量一致。字符串使用引号（单、双均可），如：yy将建立字符变量y。3、逻辑变量逻辑变量中的元素是：TRUE(或简写为T)、FALSE(F)、NA(表示缺省)。请注意必须都是大写。逻辑变量可以直接输入，如：z-c(TRUE,FALSE,F,T,NA)可建立逻辑变量z。1可以表示T，0可以表示F。所以可以直接写成：z-c(1，0，0，1,NA)在计算时，T和F也是当作1和0来用的。二、建立数据集在上面，我们建立起了三个变量：x，y，z。如果要建立起一个三变量的数据集，即含有x、y、z的二维表，则可使用data.frame()函数。当然，在这种情况下，要求变量的长度是一样的。ddtest便把数据集读入，并命名为test。注：文件的后缀不必一定要.txt，关键文件要为纯文本，里面不能有特殊格式符。如果数据文件中没有第一列记录序号，如：Price Floor Area Rooms Age Cent.heat52.00 111.0 830 5 6.2 no54.75 128.0 710 5 7.5 no57.50 101.0 1000 5 4.2 no57.50 131.0 690 6 8.8 no59.75 93.0 900 5 1.9 yes则命令语句为：read.table(“c:/test.txt”, header=TRUE)-test系统根据每个变量第一个值的类型，自动识别变量类型，如以上数据集中，除最后一个“Cent.hear”是字符变量，其他均为数值变量。2、 读入其他格式的数据库要读入其他格式数据库，必须先安装“foreign”模块。它不属于R的8个内在模块，需要在使用前安装。安装方法很简便，只需键入命令：library(foreign)即可。Epi info: R可读入epi5和6的数据库。命令如下： read.epiinfo(“文件名.rec”) 即可。 要给数据集一个名字，则是： read.epiinfo(“文件名.rec”)-名称 例如： read.epiinfo(“d:/ttt.rec”)-tttStata: R可读入Stata5,6,7的数据库。命令如下： read.dta(“文件名.dta”)SAS： 对于SAS，R只能读入SAS Transport format(XPORT)文件。所以需要把普通的SAS数据文件（.ssd和.sas7bdat）转换为Transport format(XPORT)文件,再用命令：read.xport( )SPSS：read.spss( )可读入SPSS数据文件。读入数据文件后，使用 数据集名$变量名 ，可使用各个变量。如：mean(ttt$age)便是计算数据集ttt中的变量age的均数。一、T检验1、 样本均数与总体均数的假设检验例如：样本如下：12，15，32，12，15，16，28，56，42总体均数为：20程序如下：c(12，15，32，12，15，16，28，56，42)-x建立变量xt.test(x,mu=20,alternative=”two.sided”, conf.level=0.95)T检验 注：以上程序中，x为数据变量；mu=20为总体均数；alternative=”two.sided”表示行双侧检验，如果alternative=”less”，为的单侧检验，反之，alternative=”greater”，为的单侧检验；conf.level=0.95指检验水准为0.05。以上程序结果如下： One Sample t-testdata: x t = 1.0312, df = 8, p-value = 0.3326alternative hypothesis: true mean is not equal to 20 95 percent confidence interval:13.40659 37.26008 sample estimates:mean of x 25.33333输出假设检验结果、备择假设、总体均数95可信区间和样本均数。2、 配对资料T检验注：以下都假设ctest模块已经加载。如资料如下：甲：11，16，35，24，16，32，16，24乙：13，21，23，16，24，15，32，16程序如下： c(11,16,35,24,16,32,16,24)-x c(13,21,23,16,24,15,32,16)-y t.test(x，y，alternative=”two.sided”,paired=TRUE, conf.level=0.95)结果如下： Paired t-testdata: x and y t = 0.4401, df = 7, p-value = 0.6731alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval:-7.652913 11.152913 sample estimates:mean of the differences 1.75输出假设检验结果、备择假设、两样本差值均数的可信区间和两样本差值的均数。3、 两样本成组设计的T检验3.1 如两样本资料如下： 甲：15，12，15，13，22，26，12，14 乙：23，15，24，23，16，25 程序如下： c(15,12,15,13,22,26,12,14)-x c(23,15,24,23,16,25)-y var.test(x,y)方差齐性检验 t.test(x,y,alternative=”two.sided”,paired=FALSE,var.equal=TRUE,conf.level=0.95) 结果如下： F test to compare two variancesdata: x and y F = 1.3896, num df = 7, denom df = 5, p-value = 0.7414alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval:0.2027743 7.3445113 sample estimates:ratio of variances 1.389628方差齐性检验结果为两样本方差齐性，所以在下面T检验中var.equal=TURE。 Two Sample t-testdata: x and y t = -1.8792, df = 12, p-value = 0.08471alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval:-10.5271616 0.7771616 sample estimates:mean of x mean of y 16.125 21.000T检验结果相信大家都能看懂了。3.2 如果数据是如下输入的：x 15，12，15，13，22，26，12，14，23，15，24，23，16，25group 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2则程序如下：c(15，12，15，13，22，26，12，14，23，15，24，23，16，25)-xc(1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2)-groupt.test(xgroup=1,xgroup=2,)二、秩和检验(使用wilcox.test方法)1、 单变量秩和检验样本资料如下：15，12，14，32，12，15，11。已知总体均数为20。程序如下： c(15,12,14,32,12,15,11)-x wilcox.test(x,mu=20,alternative=”two.sided”, conf.level = 0.95)结果如下：Wilcoxon signed rank test with continuity correctiondata: x V = 7, p-value = 0.2702alternative hypothesis: true mu is not equal to 20输出p值为0.2702。 2、 配对设计资料的秩和检验资料如下：甲：12,14,15,12,21,31,26,21 乙：21,32,15,21,12,14,12,15程序如下： c(12,14,15,12,21,31,26,21)-x c(21,32,15,21,12,14,12,15)-y wilcox.test(x,y,paired=TRUE, conf.level = 0.95)结果如下： Wilcoxon signed rank test with continuity correctiondata: x and y V = 15, p-value = 0.9322alternative hypothesis: true mu is not equal to 03、 成组设计资料的秩和检验如资料同上，非配对设计：甲：12,14,15,12,21,31,26,21 乙：21,32,15,21,12,14,12,15程序如下： c(12,14,15,12,21,31,26,21)-x c(21,32,15,21,12,14,12,15)-y wilcox.test(x,y,paired=FALSE, conf.level = 0.95)结果如下： Wilcoxon rank sum test with continuity correctiondata: x and y W = 34.5, p-value = 0.8305alternative hypothesis: true mu is not equal to 0T检验、秩和检验小结：T检验、秩和检验的参数格式大致是一样的。t.test ,wilcox.test(x,(y), alternative = (two.sided, less, greater),mu = , paired = (FALSE,TRUE), conf.level = 0.95) 可以只有一个变量x，也可以有两个，alternative指定双侧还是单侧检验，mu指定总体均数，paired指定是否配对，conf.level指定检验水准。三、方差分析1、 完全随机设计资料的方差分析如资料如下：A 12,15,21,31,21,15B 21,22,15,15,21C 24,15,12,15,32,15,12程序如下： c(12,15,21,31,21,15,21,22,15,15,21,24,15,12,15,32,15,12)-x factor(c(rep(1,6),rep(2,5),rep(3,7)-group建立分组变量group anova(lm(xgroup)方差分析注：rep(1,6)函数表示：将1重复6次。factor()函数则是表明group是一个分组变量lm(xgroup)表示建立x为因变量，group为自变量的线性模型 anova(lm(xgroup)表示将线性模型中的方差表输出。 结果如下： Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(F)gr 2 14365.5 7182.8 1.8869 0.2132Residuals 8 30453.2 3806.7 输出方差分析结果和残差结果。 2、 配伍组设计资料的方差分析资料如下：A B C Dgroup1 5.27 5.27 4.94 4.61 2 5.27 5.22 4.88 4.66 3 5.88 5.83 5.38 5.00 4 5.44 5.38 5.27 5.005 5.66 5.44 5.38 4.88 6 6.22 6.22 5.61 5.22 7 5.83 5.72 5.38 4.88 8 5.27 5.11 5.00 4.44程序如下： x-scan() 5.27 5.27 4.94 4.61 5.27 5.22 4.88 4.66 5.88 5.83 5.38 5.00 5.44 5.38 5.275.00 5.66 5.44 5.38 4.88 6.22 6.22 5.61 5.22 5.83 5.72 5.38 4.88 5.27 5.115.00 4.44g1-factor(c(rep(1:4),8)g2-factor(c(rep(1,4),rep(2,4),rep(3,4),rep(4,4),rep(5,4),rep(6,4),rep(7,4),rep(8,4)anova(lm(xg1+g2)程序解释： xF) g1 3 2.90438 0.96812 77.310 1.632e-11 *g2 7 2.49800 0.35686 28.497 2.379e-09 *Residuals 21 0.26297 0.01252 输出了g1、g2的方差分析结果。四、卡方检验（使用chisq.test方法）如资料如下： A B甲 15 12乙 23 12丙 24 21 程序如下： matrix(c(15,23,24,12,12,21),nc=2)-x 建立两维表nc指列数 chisq.test(x,correct=TRUE) correct表示是否矫正结果如下： Pearsons Chi-squared testdata: x X-squared = 1.3229, df = 2, p-value = 0.5161五、线性回归（采用lm()函数）1、如资料如下：x 820,780,720,867,690,787,934,679,639,820y 165,158,130,180,134,167,186,145,120,158 对y，x做线性相关回归分析，程序如下： c(820,780,720,867,690,787,934,679,639,820)-x c(165,158,130,180,134,167,186,145,120,158)-y plot(x,y)画x,y的散点图 lm(yx)-rr做相关回归分析，结果存在rr中 summary(rr)全部显示rr内容输出结果如下：首先是散点图：summary(rr)输出结果：Call:lm(formula = y x)Residuals: Min 1Q Median 3Q Max -12.4065 -4.2977 -0.6728 4.3013 11.6912 Coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) -17.3575 22.2644 -0.780 0.458 x 0.2219 0.0286 7.759 5.44e-05 *-Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1 Residual standard error: 7.879 on 8 degrees of freedomMultiple R-Squared: 0.8827, Adjusted R-squared: 0.868 F-statistic: 60.2 on 1 and 8 DF, p-value: 5.44e-05可见输出了残差情况、截距、x的系数及假设检验结果、决定系数（R-Squared）和方差分析结果（最后一行）。2、如果资料是已经输入在一个数据集中的，那么分析方法如下。如上资料，设已经输入在一个数据集dd中：程序如下： lm(yx,data=dd)-rr summary(rr)3、多元回归分析中的方程式表达 方程式：yx 表示x为自变量 方程式：yx+z 表示x,z为自变量，探讨的是x,z的主效应 方程式：yx+z+x:z 表示探讨x,z的主效应及其交互作用 方程式：yx*z 意义同x+z+x:z 方程式：y(x+z+m)2 表示x+z+m+x:z+x:m+z:m(n表示包含所有n阶以下的交互作用) 方程式：yx%in%z 表示x包含在z中，即等同于z+z:x 方程式；y(x+z+m)2-x:z 表示从(x+z+m)2中去掉x:z以上方法也适用于其他一些多元分析方法。作者：陶立波一、R的常用函数（一）、数据操作函数（二）、常用统计函数（三）、制图函数二、R的程序编制（一）、R程序中的控制结构（二）、自定义函数R是一个以函数为核心的统计软件。在R中，所有的功能几乎都是以函数的形式提供的。各种数据操纵、统计处理，都表现为：函数名：（参数1，参数2，）的形式。而通过编制自己的函数，便可以方便的扩展R的功能。一、R的常用函数（一）、数据操作函数1、c( )此函数用于建立向量（vector）。如：c(1,5,32,2) 建立向量(1,5,32,2)c(21,5:8) 建立向量(21,5,6,7,8,)。5:8表示5至8。c(“dte”,”eri”,sdi”) 建立字符串向量(“dte”,”eri”,sdi”)。c(.)-x 表示建立一个向量，其名称为x；也可以理解为把一个向量赋值给变量x。 2、rep( )重复函数。rep(1,5)表示把1重复5次。如：c(rep(2,4),3,2) 表示建立向量(2,2,2,2,3,2)。3、matrix( )此函数用于建立矩阵。如：matrix(c(1,2,3,5,4,6),nc=2) 将建立3*2的矩阵，其中nc=2指列数为2。 结果如下： ,1 ,21, 1 52, 2 43, 3 64、data.frame( )此函数用于建立数据集。例如：程序如下： c(12,15,14,12)-x c(“bb”,”cc”,”dd”,”ee”)-y c(0,0,1,1)-group 已经建立了3个变量：x,y,group data.frame(x,y,group)将建立如下数据集： 5、scan( )在数据较多的时候，用c( )函数不太方便。此时，可用scan( )函数。如：scan()-x 12 15 23 15 12 14 12 14 12 54 12 . .便可把下面的数据输入x变量中，数据块可以任意长度、在任意地方回车另起一行，数据之间以空格相互分隔。6、read.?( )这是一组读入其他格式数据文件的函数。常用的有：read.epiinfo（读入epi info数据文件）、read.dta（读入stata数据文件）等。如：read.epiinfo(“d:/test.rec”)-x 便读入test.rec文件，并命名为x。需要注意的是，在使用这些读入函数前，必须先加载foreign模块，因为这些数据读入函数不在缺省的模块中。加载foreign模块的命令是library(foreign)。7、attch( )和detach( )指定和撤销可用数据集。如：attch(dd)这样，在以后程序中，便可以直接应用dd中的变量。detach( )则相反。（二）、常用统计函数1、t.test( )用于对一个或两个样本做T检验。格式：t.test(x, y ,alternative = c(two.sided, less, greater), mu = ?,paired = (FALSE,TRUE), var.equal = (FALSE,TRUE), conf.level = .) 2、wilcox.test( )用于对一个或两个样本做wilcoxon秩和检验。格式：wilcox.test(x, y ,alternative = c(two.sided, less, greater), mu = ?,paired = (FALSE,TRUE), conf.level = .)3、chisq.test( )格式：chisq.test(x,correct=(FALSE,TRUE)注：x是一个矩阵，可采用matrix( )建立。格式2：chisq.test(x,y,correct=.)x,y是两个变量，对其交叉表做卡方检验。4、lm( )用于线性回归分析。如：lm(yx) lm(yx,data=dd) 如果y,x在数据集dd中。以上具体内容，请见前面章节。（三）、制图函数1、 plot(x,y)散点图。2、 hist(x)直方图。 3、 boxplot(x)箱式图。 等等。二、R的程序编制（一）、R程序中的控制结构1、 if else.结构语法：if (条件) 语句1 else 语句2如：x=2if (x1) y=1 else y=2得到：y为2。If后面的条件，可以用“&”（与）、“|”（或）连接，成为复合条件。如：if (x2) z=1 else z=2语句可以是多条语句，这时，需要用，语句之间用分号相隔。如：if (xx y=0 for(x in 1:5) y=y+1得到：y为5。3、while 循环结构语法：while(条件) 语句例：x-0y-0while(x5) y=y+2;x=x+1得到：x为5，y为10。在语句中使用“break”语句，便可中止循环。例：x-0y-0while(x6) break 将得到：x为4，y为8。在语句中使用“next”语句，将停止当前循环，直接进入下一轮循环。例：x-0y-0while(x3) next;y=y+2将得到：x为5，y为6。因为当x大于3后，y就不再累加了。（二）、自定义函数创制自己的函数方法如下：函数名-function(参数1，参数2，) 语句例如：创建一个函数ap，其功能为将参数的和乘以参数的差。 程序如下：ap a c(3.3, 3.4, 3.8) - b # Can also assign pointing right assign(c, c(a, 4.0, b) # Or explicitly to a variable name c # Concatenation flattens arguments1 3.1 4.2 2.7 4.1 4.0 3.3 3.4 3.8 1/c # Operate on each element of vector1 0.3225806 0.2380952 0.3703704 0.2439024 0.2500000 0.3030303 0.29411768 0.2631579 a * b # Cycle shorter vector b (but warn)1 10.23 14.28 10.26 13.53Warning message:longer object length is not a multiple of shorter object length in: a * b a+1 # 1 is treated as vector length 11 4.1 5.2 3.7 5.1下面有关于 R 语法的注标（indexing）、分割、可选参数以及其他内容的更多示例。R shell 提示符 尤其是如果您安装了 GNU readlines 是一个极好的用于研究的界面。工作时别忘记了通过 help(function) 命令来学习更多内容（您也可以使用 ?function）。Python shell 用户将发现 R shell 实在是很熟悉 而且将会高度评价两者的实用程序。 回页首温度数据集近一年以来，Brad 一直在通过他房子里面以及附近的四个温度计来收集温度数据，并使用 GnuPlot（参阅 参考资料 以获得关于 Gnuplot 更多资料的链接）将不断变化的读数窗口自动编译为可以通过 Web 访问的图形。尽管此类 hackerish 的数据收集不能真正地用于任何主要的科学用途，但是，它拥有很多 类似于 科学数据的极好特性。 每三分钟收集一次数据，这样，一年中就会有非常多数据点（四个测量位置大约共有 750,000 个）。由于温度计、传输通路或者进行记录的计算机的各种失误，有一些数据会丢失。在较为少见的一些情况下，据了解，由于计时错误，会导致单线传输通路会发生同时读取的问题。换句话说，Brad 的温度数据看起来与 相当好 的实际科学数据很类似，但是还是有一些小问题，也并不是理想的。 读取数据 温度数据收集到四个以收集位置命名的独立数据文件中，每个的格式如下： 清单 2. 初始温度数据文件的格式2003 07 25 16 04 27.5000002003 07 25 16 07 27.30000020