R语言学习

#层次聚类Data=iris[,-5];Means=sapply(data,mean);SD=sapply(data,sd); dataScale=scale(data,center=means,scale=SD);Dist=dist(dataScale,method=”euclidean”);heatmap(as.matrix(Dist),labRow=FALSE,labCol=FALSE); clusteModel=hclust(Dist,method=”ward”);result=cutree(clusteModel,k=3);table(iris[,5],result);plot(clusteModel);library(fastcluster); # kuaisu cengcijuleiclusteModel=hclust(Dist,method=”ward”);library(proxy);res=dist(data,method=”cosine”);x=c(0,0,1,1,1,1);y=c(1,0,1,1,0,1);dist(rbind(x,y),method=”Jaccard”);x=c(0,0,1.2,1,0.5,1,NA);y=c(1,0,2.3,1,0.9,1,1);d=abs(x-y);Dist=sum(d[!is.na(d)])/6;# k-means聚类clusteModel=kmeans(dataScale,centers=3,nstart=10); class(clusteModel);library(proxy);library(cluster);clustModel=pam(dataScale,k=3,metric=”Mahalanobis”); clustModel$medoidstable(iris$Species,clustModel$clustering);par(mfcol=c(1,2));plot(clustModel,which.plots=2,main=””);Plot(clustModel,which.plots=1,main=””);library(devtools);install_github(“lijian13/rinds”);rinds::bestCluster(dataScale,2:6);library(fpc);pka=kmeansruns(iris[,1:4],krange=2:6,critout=TRUE,runs=2,criterion=”asw”);#基于密度的聚类x1=seq(0,pi,lenth.out=100);y1=sin(x1)+0.1*rnorm(100);x2=1.5+seq(0,pi,length.out=100);y2=cos(x2)+0.1*rnorm(100);data=data.frame(c(x1,x2),c(y1,y2));names(data)=c(“x”,”y”);model1=kmeans(data,centers=2,nstart=10);library(“fpc”);model2=dbscan(data,eps=0.3,MinPts=4);#自组织映射library(kohonen);data=as.matrix(iris[,-5]);somModel=som(data,grid=somgrid(15,10,”hexagonal”));plot(somModel,ncolors=10,type=”dist.neighbours”);irisclass=as.numeric(iris[,5]);plot(somModel,type=”mapping”,labels=irisclass,col=irisclass+3,main=”mapping plot”);#主成分分析library(FactoMineR);data(decathlon);head(decathlon,n=2);pca1=princomp(decathlon[,1:10]);plot(pca1,type="line");res.pca=PCA(decathlon,quanti.sup=11:12,quali.sup=13);#对应分析library(MASS);data(caith);biplot(corresp(caith,nf=2),xlim=c(-0.6,0.8));#多元分析的可视化library(car);data(mpg,package="ggplot2");scatterplotMatrix(mpg[,c('displ',"cty","hwy")],diagonal="histogram",ellipse=TRU E);library(corrplot);data(mtcars);M=cor(mtcars);corrplot(M,order="hclust");#Logistic回归set.seed(1);b0=1;b1=2;b2=3;x1=rnorm(1000);x2=rnorm(1000);z=b0+b1*x1+b2*x2;pr=1/(1+exp(-z));y=rbinom(1000,1,pr);plotdata2=data.frame(x1,x2,y=factor(y));library(ggplot2);p2=ggplot(data=plotdata2,aes(x=x1,y=x2,color=y))+geom_point();print(p2);data=data.frame(x1,x2,y);model=glm(y~.,data=data,family="binomial");summary(model);w=model$coef;inter=-w[1]/w[3];slope=-w[2]/w[3];plotdata3=data.frame(cbind(x1,x2),y=factor(y));p3=ggplot(data=plotdata3,aes(x=x1,y=x2,color=y))+geom_point()+geom_abline(inter cept=inter,slope=slope);predict(model,newdata=list(x1=1,x2=3),type="response");snafile=system.file("examples","sna","lijian001.txt",package="rinds");snadf=read.table(snafile,header=FALSE,stringsAsFactors=FALSE);head(snadf)library(igraph);snaobj=graph.data.frame(snadf,directed=FALSE);class(snaobj)vcount(snaobj);ecount(snaobj);neighbors(snaobj,6,mode="all");degree(snaobj,v=6);betweenness(snaobj,v=6,directed=FALSE);closeness(snaobj,v=6);page.rank(snaobj,vids=6);similarity.dice(snaobj,vids=c(6,7));snaclass=munity(snaobj,steps=5);cl=snaclass$membership;V(snaobj)$color=rainbow(max(cl))[cl];V(snaobj)$bte=betweenness(snaobj,directed=FALSE);V(snaobj)$size=5;V(snaobj)[bte>=1800]$size=15;V(snaobj)$label=NA;V(snaobj)[bte>=1800]$label= V(snaobj)[bte>=1800]$name;plot(snaobj,layout=layout.fruchterman.reingold,vertex.size=V(snaobj)$size,vertex.color= V(snaobj)$color,bel= V(snaobj)$label,bel.cex=V(snaobj)$cex,edge.color=grey(0.5),edge.arrow.mode="-");用caret包对数据清洗并进行回归树预测set.seed(1)data(PimaIndiansDiabetes2,package='mlbench')data=PimaIndiansDiabetes2library(caret)library(caret)preProcV alues=preProcess(data[,-9],method=c('center','scale')) scaleddata=predict(preProcV alues,data[,-9])preProcbox=preProcess(scaleddata,method=c('YeoJohnson'))boxdata=predict(preProcbox,scaleddata)preProcimp=preProcess(boxdata,method='bagImpute')procdata=predict(preProcimp,boxdata)procdata$class=data[,9]library(rpart)rpartModel=rpart(class~.,data=procdata,control=rpart.control(cp=0))cptable=as.data.frame(rpartModel$cptable)cptable$errsd=cptable$xerror+cptable$xstdcpvalue=cptable[which.min(cptable$errsd),'CP']pruneModel=prune(rpartModel,cpvalue)library(rpart.plot)rpart.plot(pruneModel)pre=predict(pruneModel,procdata,type='class')preTable=table(pre,procdata$class)accuracy=sum(diag(preTable))/sum(preTable)write.table(iris,file='C:/Program Files/R/zhangfuchang.csv',sep=',') data=read.table(file='C:/Program Files/R/zhangfuchang.csv',sep=',')write.table(procdata,file='C:/Program Files/R/zhangfuchang.csv',sep=',') procdata=read.table(file='C:/Program Files/R/zhangfuchang.csv',sep=',')回归树回归代码rpartModel=rpart(class~.,data=procdata,control=rpart.control(cp=0.01),parms=list(loss=matrix(c( 0,5,1,0),2)))pre=predict(rpartModel,procdata,type='class')preTable=table(pre,procdata$class)accuracy=sum(diag(preTable))/sum(preTable)用分类回归数分类，并用10重交叉验证的R代码procdata=read.table(file='C:/Program Files/R/zhangfuchang.csv',sep=',')num=sample(1:10,nrow(procdata),replace=TRUE)res=array(0,dim=c(2,2,10))n=ncol(procdata)for (i in 1:10){train=procdata[num!=i,]test=procdata[num==i,]model=rpart(class~.,data=train,control=rpart.control(cp=0.1))pre=predict(model,test[,-n],type='class')res[,,i]=as.matrix(table(pre,test[,n]))}table=apply(res,MARGIN=c(1,2),sum)sum(diag(table)/sum(table))贝叶斯分类library(MASS)library(klaR)nbModel=NaiveBayes(class~.,data=procdata,usekernel=FALSE,fL=1)plot(nbModel,vars='glucose',legendplot=TRUE)plot(nbModel,vars='pressure',legendplot=TRUE)支持向量机分类library(devtools)install_github("lijian13/rinds")library(rinds)data(LMdata,package="rinds")library(kernlab)model=ksvm(y~.,data=LMdata$SVM,kernel='rbfdot',C=1)plot(model,data=LMdata$SVM)用caret包中的train函数分别进行回归树、贝叶斯分类、最近邻分类、神经网络分类、支持向量机分类、随即森林分类回归树：library(caret)library(caret)library(e1071)fitControl=trainControl(method="repeatedcv",number=10,repeats=3)tunedf=data.frame(.cp=seq(0.001,0.1,length=10))treemodel=train(x=procdata[,-9],y=procdata[,9],method='rpart',trControl=fitControl,tuneGrid=tun edf)plot(treemodel)贝叶斯：library(klaR)library(caret)fitControl=trainControl(method="repeatedcv",number=10,repeats=3)tunedf=data.frame(.fL=1,.usekernel=TRUE)nbmodel=train(x=procdata[,-9],y=procdata[,9],method='nb',trControl=fitControl,tuneGrid=tunedf) densityplot(nbmodel)最近邻：library(caret)fitControl=trainControl(method="repeatedcv",number=10,repeats=3)tunedf=data.frame(.k=seq(3,20,by=2))knnmodel=train(x=procdata[,-9],y=procdata[,9],method='knn',trControl=fitControl,tuneGrid=tune df)plot(knnmodel)神经网络：library(caret)library(nnet)fitControl=trainControl(method="repeatedcv",number=10,repeats=3)tunedf=expand.grid(.decay=0.1,.size=5:10,.bag=TRUE)nnetmodel=train(class~.,data=procdata,method='avNNet',trControl=fitControl,trace=FALSE,linou t=FALSE,tuneGrid=tunedf)plot(nnetmodel)支持向量机library(caret)library(kernlab)fitControl=trainControl(method="repeatedcv",number=10,repeats=3)tunedf=data.frame(.C=seq(0,1,length=11))svmmodel=train(class~.,data=procdata,method='svmRadialCost',trControl=fitControl,tuneGrid=tu nedf)plot(svmmodel)随即森林：library(caret)library(randomForest)fitControl=trainControl(method="repeatedcv",number=10,repeats=3)rfmodel=train(class~.,data=procdata,method='rf',trControl=fitControl,tuneLength=5)plot(rfmodel)varImpPlot(rfmodel$finalModel)partialPlot(rfmodel$finalModel,procdata[,-9],"mass",which.class="pos")importance(rfmodel$finalModel)rm(list=ls())x=rnorm(1e4)*2+6y=0.5*x+1y=y+rnorm(length(y))smoothScatter(x,y)model=lm(y~x)abline(model,lwd=2)text(1,8,"R^2=0.XXXX")R语言作单样本的Wilcoxon秩和检验：x=c(4.12,5.81,7.63,9.74,10.39,11.92,12.32,12.89,13.54,14.45) wilcox.test(x-8,alt="greater")R语言作成对样本的Wilcoxon秩和检验：algae=read.table("Analysis.txt",header=F,dec=".",s=c("season","size","speed","mxPH"," mnO2","Cl","NO3","NH4","oPO4","PO4","Chla","a1","a2","a3","a4","a5","a6","a7"),na.strings= c("XXXXXXX"))library(DMwR)data(algae)找出那些样本含有NA algae[!complete.cases(algae),]求出含有NA的样本总个数nrow(algae[!complete.cases(algae),])剔除含有NA的所有样本algae=na.omit(algae)每个样本含有多少个NA apply(algae,1,function(x) sum(is.na(x)))library(DMwR)data(algae)algae=algae[-manyNAs(algae),]clean.algae=knnImputation(algae,k=10)lm.a1=lm(a1~.,data=clean.algae[,1:12])anova(lm.a1)library(DMwR)data(algae)library(rpart)algae=algae[-manyNAs(algae),]rt.a1=rpart(a1~.,data=algae[,1:12])prettyTree(rt.a1)library(DMwR)data(algae)algae=algae[-manyNAs(algae),]clean.algae=knnImputation(algae,k=10)cv.rpart=function(form,train,test,...){m=rpartXse(form,train,...)p=predict(m,test)mse=mean((p-resp(form,test))^2)c(nmse=mse/mean((mean(resp(form,train))-resp(form,test))^2)) }cv.lm=function(form,train,test,...){m=lm(form,train,...)p=predict(m,test)p=ifelse(p<0,0,p)mse=mean((p-resp(form,test))^2)c(nmse=mse/mean((mean(resp(form,train))-resp(form,test))^2)) }res=experimentalComparison(c(dataset(a1~.,clean.algae[,1:12],"a1")),c(variants("cv.lm"),variants("cv.rpart",se=c(0,0.5,1))), cvSettings(3,10,1234))。