【另附】:R语言简明笔记系列1.1. 描述性统计分析1.1.1. 描述性统计量的计算1.1.1.1. summary()

>vars<-c("mpg","hp","wt")>head(mtcars[vars])mpghpwtMazdaRX421.01102.620MazdaRX4Wag21.01102.875Datsun71022.8932.320Hornet4Drive21.41103.215HornetSportabout18.71753.440Valiant18.11053.460>summary(mtcars[vars])mpghpwtMin.:10.40Min.:52.0Min.:1.5131stQu.:15.431stQu.:96.51stQu.:2.581Median:19.20Median:123.0Median:3.325Mean:20.09Mean:146.7Mean:3.2173rdQu.:22.803rdQu.:180.03rdQu.:3.610Max.:33.90Max.:335.0Max.:5.424</pre>>**`summary()`**函数提供最小值、最大值、四分位数和数值型变量的均值,以及因子向量和逻辑型向量的频数统计1.1.1.2. describe() 【Hmisc包】

>library(Hmisc)载入需要的程辑包:grid载入需要的程辑包:lattice载入需要的程辑包:survival载入需要的程辑包:Formula载入需要的程辑包:ggplot2载入程辑包:‘Hmisc’Thefollowingobjectsaremaskedfrom‘package:base’:format.pval,round.POSIXt,trunc.POSIXt,units>describe(mtcars[vars])mtcars[vars]3Variables32Observations------------------------------------------------------------------------------------------------mpgnmissinguniqueInfoMean.05.10.25.50.75.90.9532025120.0912.0014.3415.4319.2022.8030.0931.30lowest:10.413.314.314.715.0,highest:26.027.330.432.433.9------------------------------------------------------------------------------------------------hpnmissinguniqueInfoMean.05.10.25.50.75.90.95320221146.763.6566.0096.50123.00180.00243.50253.55lowest:5262656691,highest:215230245264335------------------------------------------------------------------------------------------------wtnmissinguniqueInfoMean.05.10.25.50.75.90.953202913.2171.7361.9562.5813.3253.6104.0485.293lowest:1.5131.6151.8351.9352.140,highest:3.8454.0705.2505.3455.424------------------------------------------------------------------------------------------------1.1.1.3. stat.desc() 【pastecs包】

>library(pastecs)载入需要的程辑包:boot>stat.desc(mtcars[vars])mpghpwtnbr.val32.000000032.000000032.0000000nbr.null0.00000000.00000000.0000000nbr.na0.00000000.00000000.0000000min10.400000052.00000001.5130000max33.9000000335.00000005.4240000range23.5000000283.00000003.9110000sum642.90000004694.0000000102.9520000median19.2000000123.00000003.3250000mean20.0906250146.68750003.2172500SE.mean1.065424012.12031730.1729685CI.mean.0.952.172946524.71955010.3527715var36.32410284700.86693550.9573790std.dev6.026948168.56286850.9784574coef.var0.29998810.46740770.30412851.1.1.4. describe() 【psych包】

>library(psych)载入程辑包:‘psych’Thefollowingobjectismaskedfrom‘package:boot’:logit>describe(mtcars[vars])varsnmeansdmediantrimmedmadminmaxrangeskewkurtosissempg13220.096.0319.2019.705.4110.4033.9023.500.61-0.371.07hp232146.6968.56123.00141.1977.1052.00335.00283.000.73-0.1412.12wt3323.220.983.333.150.771.515.423.910.42-0.020.171.1.2. 分组计算描述性统计量

在比较多组个体或观测时,关注的焦点经常是各组的描述性统计信息,而不是样本整体的描述性统计信息。同样地,在R中完成这个任务有若干种方法。我们将以获取变速箱类型各水平的描述性统计量开始。

1.1.2.1. aggregate()

>aggregate(mtcars[vars],by=list(am=mtcars$am),mean)ammpghpwt1017.11603.772124.41272.41>aggregate(mtcars[vars],by=list(am=mtcars$am),sd)ammpghpwt103.8353.90.777216.1784.10.617</pre>>由上面的分析结果可看出,am有两个值,根据am的两个值将`mtcars`数据集分为两组,得出上面的`mpg`,`hp`,`wt`的平均值以及标准差。<br>>其中,**`list(am=mtcars$am)`**的使用,将`am`列标注为一个更有帮助的列标签,而非`Group.1`。<br>>遗憾的是,aggregate()仅允许在每次调用中使用平均数、标准差这样的单返回值函数。它无法一次返回若干个统计量。要完成这项任务,可以使用**`by()`**函数。1.1.2.2. by()

>by(mtcars[vars],mtcars$am,summary)mtcars$am:0mpghpwtMin.:10.4Min.:62Min.:2.461stQu.:14.91stQu.:1161stQu.:3.44Median:17.3Median:175Median:3.52Mean:17.1Mean:160Mean:3.773rdQu.:19.23rdQu.:1923rdQu.:3.84Max.:24.4Max.:245Max.:5.42---------------------------------------------------------------------------------------------------mtcars$am:1mpghpwtMin.:15.0Min.:52Min.:1.511stQu.:21.01stQu.:661stQu.:1.94Median:22.8Median:109Median:2.32Mean:24.4Mean:127Mean:2.413rdQu.:30.43rdQu.:1133rdQu.:2.78Max.:33.9Max.:335Max.:3.571.1.2.3. 拓展1.1.2.4. summaryBy()【doBy包】

格式:summaryBy(formula, dataframe, FUN=function)

>**`formula`**支持格式:<br>>**var1+var2+var3+……+varN~groupvar1+groupvar2+……+groupvarN**<br>>**`~`**左边的变量是需要分析的数值型变量,右侧的变量是类别性的分组变量。<br>>**`function`**可为任何内建或用户自编的R函数。<pre>>library(doBy)载入需要的程辑包:survival载入程辑包:‘survival’Thefollowingobjectismaskedfrom‘package:boot’:aml>summaryBy(mpg+hp+wt~am,data=mtcars,FUN=mystats)ammpg.nmpg.meanmpg.stdevmpg.skewmpg.kurtosishp.nhp.meanhp.stdevhp.skewhp.kurtosiswt.nwt.meanwt.stdev101917.147373.8339660.01395038-0.803178319160.263253.90820-0.01422519-1.2096973193.7688950.7774001211324.392316.1665040.05256118-1.455352013126.846284.062321.359885860.5634635132.4110000.6169816wt.skewwt.kurtosis10.97592940.141567620.2103128-1.17373581.1.2.5. describe.by()【psych包】

1.1.3. 使用reshape包导出描述性统计量

>library(reshape)>dstats<-function(x)(c(n=length(x),mean=mean(x),sd=sd(x))+)>dfm<-melt(mtcars,measure.vars=c("mpg","hp","wt"),id.vars=c("am","cyl"))>cast(dfm,am+cyl+variable~.,dstats)amcylvariablenmeansd104mpg322.9000001.4525839204hp384.66666719.6553640304wt32.9350000.4075230406mpg419.1250001.6317169506hp4115.2500009.1787799606wt43.3887500.1162164708mpg1215.0500002.7743959808hp12194.16666733.3598379908wt124.1040830.76830691014mpg828.0750004.48385991114hp881.87500022.65541561214wt82.0422500.40934851316mpg320.5666670.75055531416hp3131.66666737.52776751516wt32.7550000.12816011618mpg215.4000000.56568541718hp2299.50000050.20458151818wt23.3700000.28284271.1.4. 结果可视化

分布特征的数值刻画的确很重要,但是这并不能代替视觉呈现。对于定量变量,我们有直方 图(6.3节)、密度图(6.4节)、箱线图(6.5节)和点图(6.6节)。它们都可以让我们洞悉那些依 赖于观察一小部分描述性统计量时忽略的细节。
目前我们考虑的函数都是为定量变量提供概述的。下一节中的函数则允许考察类别型变量的 分布。

1.2. 频数表和列联表1.2.1. 生成频数表1.2.1.1. 一维列联表

>mytable<-with(Arthritis,table(Improved))>mytableImprovedNoneSomeMarked421428>prop.table(mytable)ImprovedNoneSomeMarked0.50000000.16666670.3333333>prop.table(mytable)*100ImprovedNoneSomeMarked50.0000016.6666733.333331.2.1.2. 二维列联表

>attach(Arthritis)>mytable1<-xtabs(~Treatment+Improved,data=Arthritis)>mytable1ImprovedTreatmentNoneSomeMarkedPlacebo2977Treated13721>margin.table(mytable1,1)TreatmentPlaceboTreated4341>margin.table(mytable1,2)ImprovedNoneSomeMarked421428>prop.table(mytable,2)Errorinif(d2==0L){:需要TRUE/FALSE值的地方不可以用缺少值>prop.table(mytable1,2)ImprovedTreatmentNoneSomeMarkedPlacebo0.69047620.50000000.2500000Treated0.30952380.50000000.7500000>addmargins(mytable1)ImprovedTreatmentNoneSomeMarkedSumPlacebo297743Treated1372141Sum42142884>addmargins(prop.table(mytable1))ImprovedTreatmentNoneSomeMarkedSumPlacebo0.345238100.083333330.083333330.51190476Treated0.154761900.083333330.250000000.48809524Sum0.500000000.166666670.333333331.00000000>addmargins(prop.table(mytable1),2)ImprovedTreatmentNoneSomeMarkedSumPlacebo0.345238100.083333330.083333330.51190476Treated0.154761900.083333330.250000000.48809524

table(A,B);
xtabs(~A+B,data=mydata):~符号右方出现的为要进行交叉分类的变量,以+作为分隔;
margin.table(table(A,B) or xtabs.table(~A+B),no.ofvariables): 生成边际频数和比例
prop.table(table(A,B),no.ofvariables): 生成比例
addmargins(table(A,B) or xtabs.table(~A+B) or prop.table(mytable1), no.ofvariables): 增加第几个变量的合计,如果不加no.ofvariables则都加;

1.2.1.3. 三维列联表

>mytable2<-xtabs(~Treatment+Sex+Improved,data=Arthritis)>mytable2,,Improved=NoneSexTreatmentFemaleMalePlacebo1910Treated67,,Improved=SomeSexTreatmentFemaleMalePlacebo70Treated52,,Improved=MarkedSexTreatmentFemaleMalePlacebo61Treated165>ftable(mytable2)ImprovedNoneSomeMarkedTreatmentSexPlaceboFemale1976Male1001TreatedFemale6516Male725>margin.table(mytable2,1)TreatmentPlaceboTreated4341>margin.table(mytable2,2)SexFemaleMale5925>margin.table(mytable2,2)SexFemaleMale5925>margin.table(mytable2,3)ImprovedNoneSomeMarked421428>margin.table(mytable2,c(1,2,3)),,Improved=NoneSexTreatmentFemaleMalePlacebo1910Treated67,,Improved=SomeSexTreatmentFemaleMalePlacebo70Treated52,,Improved=MarkedSexTreatmentFemaleMalePlacebo61Treated165>margin.table(mytable2,c(1,3))ImprovedTreatmentNoneSomeMarkedPlacebo2977Treated13721>ftable(prop.table(mytable2,c(2,3))+)ImprovedNoneSomeMarkedTreatmentSexPlaceboFemale0.76000000.58333330.2727273Male0.58823530.00000000.1666667TreatedFemale0.24000000.41666670.7272727Male0.41176471.00000000.8333333>ftable(addmargins(prop.table(mytable2,c(1,2)),3))*100ImprovedNoneSomeMarkedSumTreatmentSexPlaceboFemale59.37500021.87500018.750000100.000000Male90.9090910.0000009.090909100.000000TreatedFemale22.22222218.51851959.259259100.000000Male50.00000014.28571435.714286100.0000001.2.1.4. 独立性检验1.2.1.4.1. 卡方独立性检验

>mytable3<-xtabs(~Treatment+Improved,data=Arthritis);mytable3ImprovedTreatmentNoneSomeMarkedPlacebo2977Treated13721>chisq.test(mytable3)Pearson'sChi-squaredtestdata:mytable3X-squared=13.055,df=2,p-value=0.0014631.2.1.4.2. Fisher精确检验

>mytable4<-xtabs(~Treatment+Improved,data=Arthritis);mytable4ImprovedTreatmentNoneSomeMarkedPlacebo2977Treated13721>fisher.test(mytable4)Fisher'sExactTestforCountDatadata:mytable4p-value=0.001393alternativehypothesis:two.sided1.2.1.4.3. Cochran-Mantel-Haenszel检验

>mytable5<-xtabs(~Treatment+Improved+Sex,data=Arthritis)>mantelhaen.test(mytable5)Cochran-Mantel-Haenszeltestdata:mytable5Cochran-Mantel-HaenszelM^2=14.632,df=2,p-value=0.0006647