多元统计分析及R语言建模_自定义函数: msaR.R

#***********************************************
#******  多元统计分析及R语言建模（第五版）******
#******  自定义函数: msaR.R               ******
#******  调用方式 source('msaR.R')        ******
#******  修改时间：王斌会 2020.2.1        ******
#***********************************************options(digits=4)  #设置4位为小数位数
par(mar=c(4,4,2,1),cex=0.8) #设置图形边际和字体大小#将数据框第一列设置为数据框行名
msa.X<-function(df){ X=df[,-1]; rownames(X)=df[,1]; X
}msa.andrews<-function(x){# x is a matrix or data frame of dataif (is.data.frame(x)==TRUE)x<-as.matrix(x)t<-seq(-pi, pi, pi/30)m<-nrow(x); n<-ncol(x)f<-array(0, c(m,length(t)))for(i in 1:m){f[i,]<-x[i,1]/sqrt(2)for( j in 2:n){if (j%%2==0) f[i,]<-f[i,]+x[i,j]*sin(j/2*t)elsef[i,]<-f[i,]+x[i,j]*cos(j%/%2*t)} }#plot(c(-pi,pi), c(min(f),max(f)), type="n", xlab="t", ylab="f(t)")plot(c(-pi,pi), c(min(f),max(f)), type="n", xlab="", ylab="")for(i in 1:m) lines(t, f[i,] , col=i)legend(2,max(f),rownames(x),col=1:nrow(x),lty=1:nrow(x),bty='n',cex=0.8)
}msa.coef.sd<-function(fm){  #标准回归系数b=fm$coef;bsi=apply(fm$model,2,sd);sibs=b[-1]*si[-1]/si[1]bs
}msa.cor.test<-function(X,diag=TRUE){p=ncol(X);if(diag){tp=matrix(1,p,p);for(i in 1:p){for(j in 1:i) tp[i,j]=cor.test(X[,i],X[,j])$stat;for(j in i:p) tp[i,j]=cor.test(X[,i],X[,j])$p.value;}cat("corr test: \n"); tp=round(matrix(tp,p,dimnames=list(names(X),names(X))),4)print(tp)#return(tp)cat("lower is t value, upper is p value \n")} else {cat("\n corr test: t value, p value \n"); if(is.matrix(X)) var=1:pelse var=names(X);for(i in 1:(p-1)){for(j in (i+1):p) cat(' ',var[i],'-',var[j],cor.test(X[,i],X[,j])$stat,cor.test(X[,i],X[,j])$p.value,"\n")}}
} msa.pca<-function(X,cor=FALSE,m=2,scores=TRUE,ranks=TRUE,sign=TRUE,plot=TRUE){  if(m<1) returnPC=princomp(X,cor=cor)Vi=PC$sdev^2Vari=data.frame('Variance'=Vi[1:m],'Proportion'=(Vi/sum(Vi))[1:m],'Cumulative'=(cumsum(Vi)/sum(Vi))[1:m])cat("\n")Loadi=as.matrix(PC$loadings[,1:m])Compi=as.matrix(PC$scores[,1:m])if(sign) for (i in 1:m)  if(sum(Loadi[,i])<0){Loadi[,i] = -Loadi[,i] Compi[,i] = -Compi[,i] }pca<-NULLpca$vars=Variif(m<=1) pca$loadings = data.frame(Comp1=Loadi) else pca$loadings = Loadi;if(scores & !ranks) pca$scores=round(Compi,4)if(scores & plot){plot(Compi);abline(h=0,v=0,lty=3)text(Compi,row.names(X)) #  par(mar=c(4,4,2,3))#   biplot(Compi,Loadi); abline(h=0,v=0,lty=3)#  par(mar=c(4,4,1,1))}if(scores & ranks){pca$scores=round(Compi,4)Wi=Vi[1:m];WiComp=Compi%*%Wi/sum(Wi)Rank=rank(-Comp)pca$ranks=data.frame(Comp=round(Comp,4),Rank=Rank)}pca
}msa.fa<-function(X,m=2,scores=TRUE,rotation="varimax",common=TRUE,ranks=TRUE){  if(m<1) returncat("\n")S=cor(X); p<-nrow(S); diag_S<-diag(S); sum_rank<-sum(diag_S)rowname = names(X)colname<-paste("Factor", 1:p, sep="")A<-matrix(0, nrow=p, ncol=p, dimnames=list(rowname, colname))eig<-eigen(S)for (i in 1:p)A[,i]<-sqrt(eig$values[i])*eig$vectors[,i]for (i in 1:p) { if(sum(A[,i])<0) A[,i] = -A[,i] }h<-diag(A%*%t(A))rowname<-c("Variance","Proportion","Cumulative")B<-matrix(0, nrow=3, ncol=p, dimnames=list(rowname, colname))for (i in 1:p){B[1,i]<-sum(A[,i]^2)B[2,i]<-B[1,i]/sum_rankB[3,i]<-sum(B[1,1:i])/sum_rank}W=B[2,1:m]*100; Vars=data.frame('Variance'=B[1,],'Proportion'=B[2,]*100,'Cumulative'=B[3,]*100)A=A[,1:m]if(rotation == "varimax" & m>1){   #cat("\n Factor Analysis for Princomp in Varimax: \n\n");VA=varimax(A); A=VA$loadings; s2=apply(A^2,2,sum); k=rank(-s2); s2=s2[k]; W=s2/sum(B[1,])*100; Vars=data.frame('Variance'=s2,'Proportion'=W,'Cumulative'=cumsum(W))rownames(Vars) <- paste("Factor", 1:m, sep="")A=A[,k]for (i in 1:m) { if(sum(A[,i])<0) A[,i] = -A[,i] }A=A[,1:m]; colnames(A) <- paste("Factor", 1:m, sep="")}fit<-NULLfit$vars<-round(Vars[1:m,],3)if(m<=1) fit$loadings <- data.frame("Factor1"=round(A,4))else fit$loadings <- round(A,4)if(common){   fit$common <- round(apply(A^2,1,sum),4)} Z=as.matrix(scale(X));PCs=Z%*%solve(S)%*%Afit$scores <- round(PCs,4)if(ranks){W=apply(fit$loadings^2,2,sum)Wi=W/sum(W);F=PCs%*%Wi fit$ranks=data.frame(Factor=round(F,4),Rank=rank(-F))} fit
}msa.KMO<-function(r){cl <- match.call()if (nrow(r) > ncol(r)) r <- cor(r, use = "pairwise")Q <- try(solve(r))if (class(Q) == as.character("try-error")) {message("matrix is not invertible, image not found")Q <- r}S2 <- diag(1/diag(Q))IC <- S2 %*% Q %*% S2Q <- Image <- cov2cor(Q)diag(Q) <- 0diag(r) <- 0sumQ2 <- sum(Q^2)sumr2 <- sum(r^2)MSAi <- colSums(r^2)/(colSums(r^2) + colSums(Q^2))kmo <- sumr2/(sumr2 + sumQ2)ans <- list(MSAi = MSAi, KMO = kmo,result = test)return(ans)
}msa.bartlett<-function(R, n = NULL, diag = TRUE){if (dim(R)[1] != dim(R)[2]) {n <- dim(R)[1]message("R was not square, finding R from data")R <- cor(R, use = "pairwise")}p <- dim(R)[2]if (!is.matrix(R)) R <- as.matrix(R)if (is.null(n)) {n <- 100warning("n not specified, 100 used")}if (diag) diag(R) <- 1detR <- det(R)#蠂2=-[n-(2p+11)/6]ln|R|;   df=p(p-1)/2statistic <- -log(detR) * (n - 1 - (2 * p + 5)/6)df <- p * (p - 1)/2pval <- pchisq(statistic, df, lower.tail = FALSE)bartlett <- list(chisq = statistic, df = df, p.value = pval)return(bartlett)
}msa.cor<-function (X, diag = TRUE){options(digits = 4)p = ncol(X)if (diag) {tp = matrix(1, p, p)for (i in 1:p) {for (j in 1:i) tp[i, j] = cor.test(X[, i], X[, j])$statfor (j in i:p) tp[i, j] = cor.test(X[, i], X[, j])$p.value}cat("corr test: \n")tp = round(matrix(tp, p, dimnames = list(names(X), names(X))), 4)print(tp)cat("lower is t value, upper is p value \n")}else {cat("\n corr test: t value, p value \n")if (is.matrix(X)) var = 1:pelse var = names(X)for (i in 1:(p - 1)) {for (j in (i + 1):p) cat(" ", var[i], "-", var[j], cor.test(X[,i],X[,j])$stat,cor.test(X[,i],X[,j])$p.value, "\n")}}
}msa.cancor<-function (x, y, pq=min(ncol(x),ncol(y)), plot = FALSE){x = scale(x)y = scale(y)n = nrow(x)p = ncol(x)q = ncol(y)ca = cancor(x, y)#cat("\n");  print(ca)r = ca$corm <- length(r)Q <- rep(0, m)P = rep(0, m)lambda <- 1for (k in m:1) {lambda <- lambda * (1 - r[k]^2)Q[k] <- -log(lambda)}s <- 0i <- mfor (k in 1:m) {Q[k] <- (n - k + 1 - 1/2 * (p + q + 3) + s) * Q[k]P[k] <- 1 - pchisq(Q[k], (p - k + 1) * (q - k + 1))}#cat("\n cancor test: \n")#print(round(data.frame(r, Q, P),4))cr=round(data.frame(CR=r, Q, P),4)cat("\n")u=as.data.frame(ca$xcoef[,1:pq]); colnames(u)=paste('u',1:pq,sep='')#print(round(u,4))v=as.data.frame(ca$ycoef[,1:pq]); colnames(v)=paste('v',1:pq,sep='')#print(round(v,4))if (plot) {u1 = as.matrix(x) %*% u[,1]v1 = as.matrix(y) %*% v[,1]plot(u1, v1, xlab = "u1", ylab = "v1")abline(lm(u1 ~ v1))}list(cor=cr,xcoef=t(round(u,4)),ycoef=t(round(v,4)))
}msa.AHP<-function(B){A=matrix(B,nrow=sqrt(length(B)),ncol=sqrt(length(B)),byrow=TRUE)print(A)m=ncol(A)ai=apply(A,1,prod)^(1/m)W=ai/sum(ai); if(m>2){AW=A%*%WL_max=sum(AW/W)/m; CI=(L_max-m)/(m-1); RI=c(0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45,1.49,1.51)CR=CI/RI[m]cat('\n    L_max=',L_max,'\n')cat('    CI=',CI,'\n')cat('    CR=',CR,'\n')if(CR<=0.1) cat('  Consistency test is OK!\n\n')else cat('    Please adjust the judgment matrix!\n')}return(W)
}

多元统计分析及R语言建模_自定义函数: msaR.R相关推荐

客户行为模型 r语言建模_客户行为建模：汇总统计的问题
客户行为模型 r语言建模 As a Data Scientist, I spend quite a bit of time thinking about Customer Lifetime Value ...
mcem r语言代码_生态学数据处理常用R语言代码
使用R来处理生态学数据越来越受到科研工作者的青睐,语义编程风格.漂亮的出图效果,能直接俘获众多用户.本文将生态学数据处理中经常会使用到的功能做个搜集整理. 本文假设读者有一些R的基础知识,对于R的编程 ...
多元统计分析及R语言建模（第五版）——第3章多元数据的直观表示课后习题
第3章多元数据的直观表示本文用到的数据可以去这个网址下下载多元统计分析及R语言建模(第5版)数据练习题 2)表3-2是2004年广东省各市高新技术产品情况.试对资料按照本章介绍的多元图示方法做直观 ...
多元统计分析及R语言建模
目录一.数据矩阵数据框及R表示 1.创建向量和矩阵 1)创建一个向量 2)创建一个矩阵 2.矩阵其他运算 1)矩阵的转置,加法减,矩阵相乘,求矩阵C的逆 2)获得矩阵对角线元素 ,创建三阶单位矩阵 ...
多元统计分析及R语言建模（第五版）——第6章判别分析课后习题
第6章判别分析文章会用到的数据请在这个网址下下载多元统计分析及R语言建模(第五版)数据练习题 1)考虑两个数据集x1 = [3 7 2 4 4 7],x2 = [6 9 5 7 4 8] (1) ...
R语言编写自定义函数计算R方、使用自助法Bootstrapping估计多元回归模型的R方的置信区间、可视化获得的boot对象、估计单个统计量的置信区间、分别使用分位数法和BCa法
R语言编写自定义函数计算R方.使用自助法Bootstrapping估计多元回归模型的R方的置信区间.可视化获得的boot对象.估计单个统计量的置信区间.分别使用分位数法和BCa法(Bootstrapp ...
R语言ggpubr包ggsummarystats函数可视化分组条形图(自定义分组颜色、添加抖动数据点jitter、误差条)并在X轴标签下方添加分组对应的统计值（样本数N、中位数、四分位数的间距iqr)
R语言ggpubr包ggsummarystats函数可视化分组条形图(自定义分组颜色.添加抖动数据点jitter.误差条error bar)并在X轴标签下方添加分组对应的统计值(样本数N.中位数med ...
R语言使用lm构建线性回归模型、并将目标变量对数化（log10）实战：可视化模型预测输出与实际值对比图、可视化模型的残差、模型预测中系统误差的一个例子、自定义函数计算R方指标和均方根误差RMSE
R语言使用lm构建线性回归模型.并将目标变量对数化(log10)实战:可视化模型预测输出与实际值对比图.可视化模型的残差.模型预测中系统误差的一个例子 .自定义函数计算R方指标和均方根误差RMSE 目 ...
R语言ggplot2可视化自定义多个图例（legend）标签之间的距离实战（例如，改变数据点颜色和数据点大小图例之间的距离）
R语言ggplot2可视化自定义多个图例(legend)标签之间的距离实战(例如,改变数据点颜色和数据点大小图例之间的距离) 目录

多元统计分析及R语言建模_自定义函数: msaR.R

多元统计分析及R语言建模_自定义函数: msaR.R相关推荐

最新文章

热门文章