R语言学习 (2)

#层次聚类

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 cengcijulei

clusteModel=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$medoids

table(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=https://www.360docs.net/doc/6b2764075.html,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,https://www.360docs.net/doc/6b2764075.html,bel= V(snaobj)$label,https://www.360docs.net/doc/6b2764075.html,bel.cex=

V(snaobj)$cex,edge.color=grey(0.5),edge.arrow.mode="-");

用caret包对数据清洗并进行回归树预测

set.seed(1)

data(PimaIndiansDiabetes2,package='mlbench')

data=PimaIndiansDiabetes2

library(caret)

library(caret)

preProcValues=preProcess(data[,-9],method=c('center','scale')) scaleddata=predict(preProcValues,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$xstd

cpvalue=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 in1: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+6

y=0.5*x+1

y=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=".",https://www.360docs.net/doc/6b2764075.html,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))