data=read.table('clipboard',T)这的里read.table是R读取外部数据的常用命令,T表示第一行是表头信息,整个数据存在名为data的变量中。另一种更方便的导入方法是利用Rstudio的功能,在workspace菜单选择“import dataset”也是一样的。
newdata=edit(data)另一种情况就是我们可能只关注数据的一部分,例如从原数据中抽取第20到30号样本的Sepal.Width变量数据,因为Sepal.Width变量是第2个变量,所以此时键入下面的命令即可:
newdata=data[20:30,2]如果需要抽取所有数据的Sepal.Width变量,那么下面两个命令是等价的:
newdata=data[,2]第三种情况是需要对数据进行一些运算,例如需要将所有样本的Sepal.Width变量都放大10倍,我们先将原数据进行一个复制,再用$符号来提取运算对象即可:
newdata=data$Sepal.Width
newdata=data3 描述统计
newdata$Sepal.Width=newdata$Sepal.Width*10
mean(data$Sepal.Width)计算分类数据Species变量的频数表和条形图
sd(data$Sepal.Width)
table(data$Species)对于一元数值数据,绘制直方图和箱线图观察其分布是常用的方法:
barplot(table(data$Species))
hist(data$Sepal.Width)对于二元数值数据,则可以通过散点图来观察规律
boxplot(data$Sepal.Width)
plot(data$Sepal.Width,Sepal.Length)如果需要保存绘图结果,建议使用Rstudio中的plot菜单命令,选择save plot as image
data=iris[,c(4,5)]下一步我们想计算不同种类花瓣的平均宽度,可以使用tapply函数,在计算前先用attach命令将data这个数据框解包以方便直接操作其变量,而不需再用$符号。
attach(data)结果如下
tapply(X=Petal.Width,INDEX=Species,FUN=mean)
setosa versicolor virginica和tapply类似的还有sapply函数,在进一步讲解前初学者还需搞清楚两种数据表现方式,即stack(堆叠数据)和unstack(非堆叠数据),上面的data就是一个堆叠数据,每一行表示一个样本。而非堆叠数据可以根据unstack函数转换而来
0.246 1.326 2.026
data.unstack=unstack(data)你应该明白这二者之间的区别了,如果要对非堆叠数据计算不同种类花瓣的平均宽度,可以利用如下函数。
head(data.unstack)
sapply(data.unstack,FUN=mean)结果是一样的,也就是说tapply对应于stack数据,而sapply对应于unstack数据
x=1:10有放回抽取则是
sample(x,size=5)
sample(x,size=5,replace=T)
x=rnorm(50,mean=10,sd=5) #随机生成50个均值为10,标准差为5的随机数为作为研究对象也可以直接利用R语言内置函数t.test
mean(x)-qt(0.975,49)*sd(x)/sqrt(50) #根据统计学区间估计公式,得到95%置信度下的区间下界
mean(x) qt(0.975,49)*sd(x)/sqrt(50) #95%置信度下的区间上界
t.test(x,conf.level=0.95)从如下结果可得95%置信区间为(9.56,12.36)
One Sample t-test2 对总体均值进行假设检验
data: x
t = 15.7301, df = 49, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
9.563346 12.364729
sample estimates:
mean of x
10.96404
t.test(x,mu=10,alternative='two.sided') #这里的原假设是总体均值(mu)为10,使用双侧检验,得到P值为0.17,可见P值不够小,不能够拒绝原假设。T检验是极为常用的检验方法,除了单样本推断之外,t.test命令还可以实现两样本推断和配对样本推断。如果要对总体比率或总体方差进行推断,可以使用prop.test和var.test。
shapiro.test(x)结果所下:
Shapiro-Wilk normality test该检验的原假设是服从正态分布,由P值为0.82可判断不能拒绝总体服从正态的假设
data: x
W = 0.9863, p-value = 0.8265
wilcox.test(x,conf.int=T) #指定conf.int让函数返回中位数的置信区间5 独立性检验(联列表检验)
wilcox.test(x,mu=1) #指定mu让函数返回中位数为10的检验结果
data(CO2) #读入内置的数据包,其中Type和Treatmen是其中两个分类变量。
chisq.test(table(CO2$Type,CO2$Treatment)) #使用卡方检验函数来检验这两个因子之间是否独立
plot(Volume~Girth,data=trees,pch=16,col='red')首先绘制了两变量的散点图,然后用lm函数建立线性回归模型,并将回归直线加在原图上,最后用summary将模型结果进行了展示,从变量P值和F统计量可得回归模型是显著的。但截距项不应该为负数,所以也可以用下面方法将截距强制为0。
model=lm(Volume~Girth,data=trees)
abline(model,lty=2)
summary(model)
model2=lm(Volume~Girth-1,data=trees)2 模型诊断
par(mfrow=c(2,2))这里左上图是残差对拟合值作图,整体呈现出一种先下降后下升的模式,显示残差中可能还存在未提炼出来的影响因素。右上图残差QQ图,用以观察残差是否符合正态分布。左下图是标准化残差对拟合值,用于判断模型残差是否等方差。右下图是标准化残差对杠杆值,虚线表示的cooks距离等高线。我们发现31号样本有较大的影响。
plot(model)
par(mfrow=c(1,1))
plot(sqrt(Volume)~Girth,data=trees,pch=16,col='red')
model2=lm(sqrt(Volume)~Girth,data=trees)
abline(model2,lty=2)
summary(model2)
plot(sqrt(Volume)~Girth,data=trees,pch=16,col='red')
model2=lm(sqrt(Volume)~Girth,data=trees)
data.pre=data.frame(predict(model2,interval='prediction'))
lines(data.pre$lwr~trees$Girth,col='blue',lty=2)
lines(data.pre$upr~trees$Girth,col='blue',lty=2)
library(DAAG)
head(anesthetic)
cdplot(factor(nomove)~conc,data=anesthetic,main='条件密度图',ylab='病人移动',xlab='麻醉剂量')
anes1=glm(nomove~conc,family=binomial(link='logit'),data=anesthetic)上面的方法是使用原始的0-1数据进行建模,即每一行数据均表示一个个体,另一种是使用汇总数据进行建模,先将原始数据按下面步骤进行汇总
summary(anes1)
anestot=aggregate(anesthetic[,c('move','nomove')],by=list(conc=anesthetic$conc),FUN=sum)得到汇总数据anestot如下所示
anestot$conc=as.numeric(as.character(anestot$conc))
anestot$total=apply(anestot[,c('move','nomove')],1,sum)
anestot$prop=anestot$nomove/anestot$total
conc move nomove total prop对于汇总数据,有两种方法可以得到同样的结果,一种是将两种结果的向量合并做为因变量,如anes2模型。另一种是将比率做为因变量,总量做为权重进行建模,如anes3模型。这两种建模结果是一样的。
1 0.8 6 1 7 0.1428571
2 1.0 4 1 5 0.2000000
3 1.2 2 4 6 0.6666667
4 1.4 2 4 6 0.6666667
5 1.6 0 4 4 1.0000000
6 2.5 0 2 2 1.0000000
anes2=glm(cbind(nomove,move)~conc,family=binomial(link='logit'),data=anestot)根据logistic模型,我们可以使用predict函数来预测结果,下面根据上述模型来绘图
anes3=glm(prop~conc,family=binomial(link='logit'),weights=total,data=anestot)
x=seq(from=0,to=3,length.out=30)
y=predict(anes1,data.frame(conc=x),type='response')
plot(prop~conc,pch=16,col='red',data=anestot,xlim=c(0.5,3),main='Logistic回归曲线图',ylab='病人静止概率',xlab='麻醉剂量')
lines(y~x,lty=2,col='blue')
联系客服