R语言与属性数据分析

1. 属性数据的描述性分析

1）属性数据形式的转化

* 将数据框中的分类数据（统计好的）转换成频数表
	# input data
	Reason<-c("操作","设备","工具","工艺","材料")
	Count<-c(22,526,292,89,47)
	data0<-data.frame(Reason,Count);data0
	# change to the table
	ftable<-table(rep(Reason,Count));data1
	names(Count)<-Reason;Count

* 将数据框中的分类数据（未统计好的）转换成频数表
	ftable<-table(data$classvar);ftable

* 频数表转换成频率表
	ftable/sum(ftable)
	prop.table(ftable)

* 列联表和概率列联表
	ctable<-table(data);ctable
	prop.table(ctable)

* 将频率表转化为数据集
	Indicator<-c(1:5) // c('A','B','C')
	data<-rep(Indicator,ftable);data

2）属性数据的描述性统计图

	barplot(Count,main='条形图')
	pie(Count,radius = 1,main='饼图')
	library(qcc)
	pareto.chart(Count, ylab = "频数",ylab2 = "累计百分比(%)",main='帕累托图')

3）属性数据的基本统计量

* 描述中心位置的统计量
	median()
	quantile( , probs = c(0.25,0.75))
	getmode()
	# the mode function
	getmode<-function(x){
   
	  tab<-table(x)
	  tab.max<-max(tab)
	  if (all(tab == tab.max)) {
   mod<-NA}
	  else if(is.numeric(x)) {
   mod<-as.numeric(names(tab)[tab == tab.max])}
	  else {
   mod<-names(tab)[tab == tab.max]}
	  names(mod)<-'mode'
	  return(mod)
	}

* 描述离散程度的统计量(离异比率、GS指数、熵)
	var.ratio<-function(x){
   
	  tab<-table(x)
	  mod<-getmode(x)
	  ratio=1-tab[names(tab)==mod]/sum(tab)
	  names(ratio)<-'variation ratio'
	  return(ratio)
	}
	GS.index<-function(x){
   
	  tab<-table(x)
	  p<-prop.table(tab)
	  index<-1-sum(p^2)
	  names(index)<-'GS index'
	  return(index)
	}
	entropy<-function(x){
   
	  #x is the raw data
	  tab<-table(x)
	  p<-prop.table(tab)
	  entropy<- (-1)*sum(p*log(p))
	  names(entropy)<-'entropy'
	  return(entropy)
	}

2. 单一属性分类数据

1）分类数据的概率检验

* 不含参数的Chi2检验
	chisq.test(table, correct = TURE/FALSE, p = p0, rescale.p = TRUE/FALSE)
* 不含参数的LRT检验
	library(RVAideMemoire)
	G.test(table,p=p0)
# 等概率检验时，p0 = rep(1/length(x2), length(x2))
---------------------------------------------------------------

* 参数的极大似然估计(MLE) (用数值方法解)
* 含一个参数
	fn_likelihood<-function(p){
   
	  ... #log-likelihood function
	}
	result<-optimize(f=fn_likelihood, interval=c(0,1), maximum = TRUE);result
	# optimize函数得到maximum和objective两个值
	p_hat<-result$maximum  # the MLE for p
	p<-c(...)
	
* 含多个参数
	fn_likelihood<-function(x){
   
	  p<-x[1]
	  q<-x[2]
	  ... #log-likelihood function
	}
	result<-optim(par = c(0.1,0.1),fn_likelihood,control = list(fnscale = -1));result
	p_hat<-result$par[1]  # the MLE for p
	q_hat<-result$par[2]  # the MLE for q
---------------------------------------------------------------

* 含参数的Chi2检验
	Chistat<-chisq.test(x, correct = FALSE, p = p, rescale.p = FALSE)$statistic
* 含参数的LRT检验
	LRTstat<-G.test(x,p = p)$statistic
---------------------------------------------------------------

* 统计决策
	# critical value
	qchisq(p=0.95, df=r-m-1,lower.tail = TRUE)
	# p value
	pchisq(Chistat,df=r-m-1,lower.tail = FALSE)
	pchisq(LRTstat,df=r-m-1,lower.tail = FALSE)

3. 四格表

1）四格表的独立性检验

* U检验（单边）
	U.stat=sqrt(n)*(n11*n22-n12*n21)/sqrt(n1_*n2_*n_1*n_2) #未修正
	U.stat=sqrt(n)*(n11*n22-n12*n21 ± n/2)/sqrt(n1_*n2_*n_1*n_2) #修正
	pvalue=pnorm(U.stat, mean = 0, sd = 1, lower.tail = FALSE)
* Chi2检验（双边）
	Chi.stat=n*(n11*n22-n12*n21)^2/(n1_*n2_*n_1*n_2) #未修正
	Chi.stat=n*(abs(n11*n22-n12*n21)-n/2)