R語言常用統(tǒng)計方法實現(xiàn)課件

常用統(tǒng)計方法用R實現(xiàn)常用統(tǒng)計方法用R實現(xiàn)R中內嵌的分布R提供了四類有關統(tǒng)計分布的函數(shù):密度函數(shù)、(累積)分布函數(shù)、分位數(shù)函數(shù)、隨機數(shù)函數(shù).它們都與分布的英文名稱(或者其縮寫)相對應.對于所給的分布名稱，加前綴“d”就是密度函數(shù)(對于離散分布,指分布律);加前綴“p”就是分布函數(shù);加前綴“q”就是分位數(shù)函數(shù);加前綴“r”就是隨機數(shù)發(fā)生函數(shù).這四類函數(shù)的第一個參數(shù)是有規(guī)律的:形為dfunc的函數(shù)為x，pfunc的函數(shù)為q,qfunc的函數(shù)為p,rfunc的函數(shù)為n(但rhyper和rwilcox是特例，他們的第一個參數(shù)為nn).若R中分布的函數(shù)名為func,則四類函數(shù)的調用格式為:1)概率密度函數(shù):dfunc(x,p1,p2,...),x為數(shù)值向量;2)(累積)分布函數(shù):pfunc(q,p1,p2,...),q為數(shù)值向量;3)分位數(shù)函數(shù):qfunc(p,p1,p2,...),p為由概率構成的向量;4)隨機數(shù)函數(shù):rfunc(n,p1,p2,...),n為生成數(shù)據的個數(shù),p1,p2,...是分布的參數(shù)值.R中內嵌的分布R提供了四類有關統(tǒng)計分布的函數(shù):密度函數(shù)、(R提供的常用分布(要加前綴)分布名稱 R名稱 選項beta beta shape1,shape2binomial binom size,probCauchy cauchy location=0,scale=1chi-sqaured(

2) chisq df,ncpexponential exp rateFisher-Snedecor(F) f df1,df2,ncpgamma gamma shape,scale=1geometric geom probhypergeometric hyper m,n,klognormal lnorm meanlog=0,sdlog=1logistic logis location=0,scale=1multinomial multinom size,probnormal norm mean=0,sd=1negativebinomial nbinom size,probPoisson pois lambdaStudent's(t) t dfuniform unif min=0,max=1Weibull weibull shape,scale=1Wilcoxon'sstatisticswilcox m,n signrank nR提供的常用分布(要加前綴)分布名稱 R名稱 選項描述性統(tǒng)計位置的度量:均值、順序統(tǒng)計量、中位數(shù)、百分位數(shù)。均值計算：若x是向量、矩陣，則mean(x)返回其全部元素均值。若要返回數(shù)組某一維的均值：apply(x,dim,mean);dim=1計算行均值，dim=2計算列均值。若x是數(shù)據框，則mean(x)返回各列的均值Mean的一般用法：

mean(x,trim=0,na.rm=FALSE)trim指定去掉x兩端數(shù)的比例；na.rm=TRUE允許有缺失值。類似有sum(x)函數(shù)可求x的和。描述性統(tǒng)計位置的度量:均值、順序統(tǒng)計量、中位數(shù)、百分位數(shù)。順序統(tǒng)計量將n個數(shù)據(觀測值)按從小到大的順序排列后,稱其為順序統(tǒng)計量.函數(shù)sort(x)給出了樣本x的順序統(tǒng)計量order()給出排序后的下標rank()給出了樣本x的秩次統(tǒng)計量x<-c(75,64,47.4,66.9,62.2,62.2,58.7,63.5)sort(x)order(x)順序統(tǒng)計量將n個數(shù)據(觀測值)按從小到大的順序排列后,稱其為中位數(shù)中位數(shù)描述數(shù)據中心位置的數(shù)字特征.大體上比中位數(shù)大或小的數(shù)據個數(shù)為整個數(shù)據的一半.對于對稱分布的數(shù)據，均值與中位數(shù)比較接近;對于偏態(tài)分布的數(shù)據，均值與中位數(shù)不同.中位數(shù)的又一顯著特點是不受異常值的影響，具有穩(wěn)健性，因此它是數(shù)據分析中相當重要的統(tǒng)計量.在R軟件中，函數(shù)median()給觀測量的中位數(shù).如x<-c(75,64,47.4,66.9,62.2,62.2,58.7,63.5)median(x)median(x,na.rm=TRUE)#若數(shù)據中有缺失值中位數(shù)中位數(shù)描述數(shù)據中心位置的數(shù)字特征.大體上比中位數(shù)大或小百分位數(shù)百分位數(shù)(percentile)是中位數(shù)的推廣.將數(shù)據按從小到大的排列后，0<p<1，它的p分位點定義為：在R軟件中,quantile()函數(shù)計算觀測量的百分位數(shù).如w<-c(75.0,64.0,47.4,66.9,62.2,62.2,58.7,63.5,66.6,64.0,57.0,69.0,56.9,50.0,72.0)quantile(w)一般用法：

quantile(x,probs=seq(0,1,0.25),na.rm=FALSE)百分位數(shù)百分位數(shù)(percentile)是中位數(shù)的推廣.將數(shù)分散程度的度量表示數(shù)據分散(或變異)程度的特征量有方差、標準差、極差、四分位極差、變異系數(shù)和標準誤等.在R軟件中，用var()和sd()計算方差、標準差：var(x,na.rm=FALSE,)sd(x,na.rm=FALSE)分散程度的度量表示數(shù)據分散(或變異)程度的特征量有方差、標準變異系數(shù)、平方和對于變異系數(shù)、校正平方和、未校正平方和等指標，需要編寫簡單的程序.變異系數(shù)CV計算:

cv<-100*sd(x)/mean(x);cv校正平方和CSS:css<-sum((x-mean(x))^2);css未校正平方和USS:uss<-sum(x^2);uss變異系數(shù)、平方和對于變異系數(shù)、校正平方和、未校正平方和等指標極差與標準誤樣本極差(記為R)的計算:R=max(x)-min(x)樣本上、下四分位數(shù)之差稱為四分位差(或半極差)，記為R1.它也是度量樣本分散性的重要數(shù)字特征，特別對于具有異常值的數(shù)據，它作為分散性具有穩(wěn)健性，因此在穩(wěn)健性數(shù)據分析中具有重要作用.半極差計算：R1=quantile(x,0.75)-quantile(x,0.25)樣本標準誤(記為sm)定義為s/sqrt(n)樣本標準誤計算：sm=sd(x)/sqrt(length(x))極差與標準誤樣本極差(記為R)的計算:R=max(分布形狀的度量偏度系數(shù)Kurtosis是刻劃數(shù)據的對稱性指標.關于均值對稱的數(shù)據其偏度系數(shù)為0.右側更分散的數(shù)據偏度系數(shù)為正,左側更分散的數(shù)據偏度系數(shù)為負.當數(shù)據的總體分布為正態(tài)分布時，峰度系數(shù)Skewness近似為0;當峰度系數(shù)為正時,兩側極端數(shù)據較多;當峰度系數(shù)為負時,兩側極端數(shù)據較少.分布形狀的度量偏度系數(shù)Kurtosis是刻劃數(shù)據的對稱性指標偏度系數(shù)Skewness樣本峰度系數(shù)sk計算程序n<-length(x)m<-mean(x)s<-sd(x)sk<-n/((n-1)*(n-2))*sum((x-m)^3)/s^3計算公式偏度系數(shù)Skewness樣本峰度系數(shù)sk計算程序峰度系數(shù)Kurtosis計算樣本峰度系數(shù)ku計算程序n<-length(x)m<-mean(x)s<-sd(x)ku<-((n*(n+1))/((n-1)*(n-2)*(n-3))*sum((x-m)^4)/s^4-(3*(n-1)^2)/((n-2)*(n-3)))計算公式峰度系數(shù)Kurtosis計算樣本峰度系數(shù)ku計算程序莖葉圖莖葉圖也是考查數(shù)據分布的重要方法R中用stem()制作莖葉圖：stem(x,scale=1,width=80,atom=1e-08)其中x是數(shù)據向量,scale控制莖葉的長度，width控制寬度,atom控制容差;如果選擇scale=2,即將10個個位數(shù)分成兩段,0?4為一段,5?9為另一段

.x<-c(25,45,50,54,55,61,64,68,72,75,75,78,79,81,83,84,

84,84,85,86,86,86,87,89,89,89,90,91,91,92,100)stem(x)stem(x,scale=2)

stem(x,scale=.5)

莖葉圖莖葉圖也是考查數(shù)據分布的重要方法直方圖直方圖是探索性數(shù)據分析的基本工具，它給出了數(shù)據的頻率分布圖形，在組距相等場合下常用寬度相等的長條矩形表示，矩形的高低表示頻率的大小在圖形上，橫坐標表示所關心變量的取值區(qū)間，縱坐標表示頻率或頻數(shù)的大小。這樣就得到頻數(shù)或頻數(shù)直方圖。圖形的形狀與我們選擇的各組區(qū)間端點有關，故選擇區(qū)間端點時我們要謹慎。R中使用函數(shù)hist()來畫直方圖，其常用的調用格式如下：

hist(x,breaks="Sturges",freq=NULL,probability=!freq,col=NULL,main=paste("Histogramof",xname),xlim=range(breaks),ylim=NULL,xlab=xname,ylab,axes=TRUE,nclass=NULL)說明:若選項breaks取向量，則用于指明直方圖區(qū)間的分割位置;若取正整數(shù),則用于指定直方圖的小區(qū)間數(shù).freq取T表示使用頻數(shù)畫直方圖,取F則使用頻率畫直方圖.probability與freq恰好相反.col用于指明小矩形的顏色.直方圖直方圖是探索性數(shù)據分析的基本工具，它給出了數(shù)據的頻率分直方圖示例已知15位學生體重ww<-c(75.0,64.0,47.4,66.9,62.2,62.2,58.7,63.5,66.6,64.0,57.0,69.0,56.9,50.0,72.0)hist(w,freq=FALSE)lines(density(w),col="blue")x<-44:76lines(x,dnorm(x,mean(w),sd(w)),col="red")執(zhí)行后繪出直方圖和密度估計曲線和正態(tài)分布的概率密度曲線直方圖示例已知15位學生體重w條形圖barplot(height,names.arg=NULL,legend.text=NULL,beside=FALSE,horiz=FALSE,col=NULL,xlab=NULL,ylab=NULL)height為表示各條形高度的向量或矩陣,names.arg為各條形名稱向量;beside控制條塊是否非堆疊式;horiz控制條形方向是否水平;col為表示各條形的顏色向量;legend.text表示圖例的文本向量pie.sales=c(0.12,0.3,0.26,0.16,0.04,0.12)pie.col=c("purple","violetred1","green3","cornsilk","cyan","white")barplot(pie.sales,col=pie.col)t=c(.20,.20,.60,.20,.30,.50,.10,.30,.6);dim(t)=c(3,3)lt=c("高中及以下","大專","本科及以上")nt=c("上海","廣州","北京")barplot(t,beside=T,xlab="城市",ylab="比例",legend.text=lt,names.arg=nt);條形圖barplot(height,names.arg=圓餅圖pie(x,labels=names(x),main=NULL,...)x表示各扇形的比例向量;labels表示各扇形的名稱向量;main為標題串。pie.sales=c(0.12,0.3,0.26,0.16,0.04,0.12)names(pie.sales)=c("Blueberry","Cherry","Apple","BostonCream","Other","VanillaCream")

pie.col=c("purple","violetred1","green3","cornsilk","cyan","white")pie(pie.sales,col=pie.col)圓餅圖pie(x,labels=names(x),頻數(shù)統(tǒng)計R中可用table(data)來進行頻數(shù)統(tǒng)計。示例1：dd=c("a","b","c","a","c","b","c","c","a");ddx=table(dd);x#頻數(shù)統(tǒng)計pie(x,labels=paste(names(x),'\n',round(100*x/sum(x),1),'%'))#繪制餅圖,paste函數(shù)可將多個字符串聯(lián)接示例2：dt1=read.table("d:/tmp/pstj1.txt")x=table(dt1);xas.data.frame(x,responseName="頻數(shù)")#轉換為表格顯示as.data.frame(x/sum(x),responseName="頻率")#顯示頻率表pie(x,labels=paste(names(x),'\n',100*x/sum(x),'%'))#圓餅圖barplot(x,names.arg=1:5,legend.text=paste(1:5,'>',names(x)))#繪制條形圖頻數(shù)統(tǒng)計R中可用table(data)來進行頻數(shù)統(tǒng)計?？ǚ綑z驗:chisq.test()chisq.test(x,y=NULL,correct=TRUE,p=rep(1/length(x),length(x)),rescale.p=F)功能：可用于列聯(lián)表檢驗和擬合優(yōu)度檢驗x是要分析的數(shù)據向量或矩陣;在做列聯(lián)表檢驗時,若x僅提供的是向量,則需要提供同樣大小的y向量;correct指示要做連續(xù)性校正;p為擬合優(yōu)度檢驗中要驗證的各分類理論頻率或頻數(shù);rescale.p指示是否要將p參數(shù)中的數(shù)據歸一化?？ǚ綑z驗:chisq.test()chisq.test(x擬合優(yōu)度檢驗示例按顏色與形狀把豌豆分為四類:黃而圓的、青而圓的、黃而有角的、青而有角的，遺傳理論告訴我們：這四類的比例為9:3:3:1。現(xiàn)觀測556個豌豆，發(fā)現(xiàn)該四類的實際數(shù)目依次為315、108、101、32，試判斷是否符合已知的遺傳規(guī)律。chisq.test(x=c(315,108,101,32),p=c(9,3,3,1),rescale.p=T)結果輸出：

Chi-squaredtestforgivenprobabilitiesdata:c(315,108,101,32)X-squared=0.47,df=3,p-value=0.9254結論：p=0.9254>0.05,接受原假設,與理論頻率無明顯差異。符合已知規(guī)律。擬合優(yōu)度檢驗示例按顏色與形狀把豌豆分為四類:黃而圓的、青而圓帶參數(shù)的分類數(shù)據的檢驗遺傳模型告知:男性正常、女性正常、男性色盲和女性色盲這四類人所占的比例分別為p/2、p2/2+p(1-p)、(1-p)/2、(1-p)2/2，其中p未知?，F(xiàn)隨機調查1000人按發(fā)現(xiàn)男性正常、女性正常、男性色盲和女性色盲各有442、514、38、6人。問調查數(shù)據是否與上述模型相符？先用最大似然法估計出參數(shù)p:求解非線性方程：f=function(p)956/p-514/(2-p)-50/(1-p)r=uniroot(f,c(0,1))(p=r$root)結果輸出：(p的值)[1]0.9129697帶參數(shù)的分類數(shù)據的檢驗遺傳模型告知:男性正常、女性正常、男續(xù)上:繼續(xù)求解的R指令#下面3條R指令計算

2值fn=1000*c(p/2,(p^2)/2+p*(1-p),(1-p)/2,((1-p)^2)/2)#理論頻數(shù)fno=c(442,514,38,6)#實際觀測頻數(shù)x2=sum((fno-fn)^2/fn);x2#卡方值3.089209#下面計算卡方檢驗概率r=length(fno);m=1;df=r-m-1#卡方維數(shù)df,觀測數(shù)r,參數(shù)m個x2p=1-pchisq(x2,df);x2p#卡方檢驗概率x2p=0.2133962結論：檢驗概率值為0.2133962>0.05,所以不拒絕原假設;我們認為關于性別和色盲的觀察數(shù)據與遺傳學的理論相符。續(xù)上:繼續(xù)求解的R指令#下面3條R指令計算2值附注：非線性方程組求解的R指令需要下載rootSolve包并裝載進內存multiroot(f,start,maxiter=100,rtol=1e-6,atol=1e-8,ctol=1e-8,useFortran=TRUE,positive=FALSE,jacfunc=NULL,jactype="fullint",verbose=FALSE,bandup=1,banddown=1,...)用法示例：library(rootSolve)model<-function(x)c(F1=x[1]^2+x[2]^2-1,F2=x[1]^2-x[2]^2+0.5)(ss<-multiroot(f=model,start=c(1,1)))結果輸出：解、函數(shù)值、迭代次數(shù)、誤差$root[1]0.50000000.8660254$f.rootF1F22.323138e-082.323308e-08$iter[1]5$estim.precis[1]2.323223e-08附注：非線性方程組求解的R指令需要下載rootSolve包并列聯(lián)表檢驗示例對63個肺癌患者和43人組成的對照組的調查數(shù)據如表,問總體中患肺癌是否與吸煙有關系?(

=0.05)compare<-matrix(c(60,32,3,11),nr=2,dimnames=list(c("cancer","normal"),c("smoke","Notsmoke")))chisq.test(compare,correct=TRUE)列聯(lián)表檢驗示例對63個肺癌患者和43人組成的對照組的調查數(shù)據交叉表(列聯(lián)表)調用格式:xtabs(formula=~.,data=parent.frame())公式的格式為:頻數(shù)變量~因子1+因子2;公式中若有點號表示數(shù)據框中除已出現(xiàn)變量外的其它變量。注意：用xtabs做卡方檢驗時，沒有連續(xù)校正的處理；而用chisq.test時默認做連續(xù)校正的處理。sex=rep(c('f','m'),c(3,3));sexoptions=rep(1:3,2);optionsnum=c(21,33,48,18,45,23);numcc=xtabs(num~sex+options);ccsummary(cc)DF<-as.data.frame(UCBAdmissions)cc=xtabs(Freq~Gender+Admit,DF);ccsummary(cc);summary(xtabs(Freq~.,DF))交叉表(列聯(lián)表)調用格式:xtabs(formula=Fisher精確檢驗fisher.test(x,y=NULL,alternative="two.sided")x是要分析的數(shù)據矩陣或向量,當x是向量時,要提供同樣大小的因子向量y;alternative指示備擇假設類型是"two.sided","greater"or"less"。數(shù)據同前例,問總體中肺癌患者吸煙的比例是否比對照組中吸煙的比例要大?(α=0.05)compare<-matrix(c(60,32,3,11),nr=2,dimnames=list(c("cancer","normal"),c("smoke","Notsmoke")))fisher.test(compare,alternative="greater")Fisher精確檢驗fisher.test(x,y=N相關分析R軟件采用用cov()函數(shù)計算協(xié)方差或協(xié)方差陣，用cor()函數(shù)計算相關矩陣(相關系數(shù))。函數(shù)cov()和cor()的使用格式為：cov(x,y=NULL,use="all.obs“,method=c("pearson","kendall","spearman"))cor(x,y=NULL,use="all.obs“,method=c("pearson","kendall","spearman"))其中x是數(shù)值型向量、矩陣或數(shù)據框.y是空值(NULL，缺省值)、向量、矩陣或數(shù)據框，但需要與x的維數(shù)相一致.

與cov和cor有關的函數(shù)還有:

cov.wt----計算加權協(xié)方差(加權協(xié)方差矩陣);cor.test---計算相關性檢驗.相關分析R軟件采用用cov()函數(shù)計算協(xié)方差或協(xié)方差陣，用c相關分析示例例為了解某種橡膠的性能，今抽取10個樣品，每個測量三項指標:硬度、變形和彈性(rubber.txt).試計算樣本均值、樣本協(xié)方差陣和樣本相關矩陣.并用Pearson相關性檢驗確認變量X1,X2,X3是否相關?rubber<-read.table("d:/rubber.txt")mean(rubber)cov(rubber)cor(rubber)cor.test(~X1+X2,data=rubber)cor.test(~X1+X3,data=rubber)cor.test(~X2+X3,data=rubber)相關分析示例例為了解某種橡膠的性能，今抽取10個樣品，每個測回歸分析案例：根據經驗，在人的身高相等的情況下，血壓的收縮壓Y與體重X1(千克)、年齡X2(歲數(shù))有關.現(xiàn)收集了13個男子的數(shù)據，見表.試建立Y關于X1，X2的線性回歸方程.估計出Y=b0+b1X1+b2X2F檢驗:

H0:b1=b2=0.T檢驗:H0:bj=0j=0,1,2回歸分析案例：根據經驗，在人的身高相等的情況下，血壓的收縮壓求解程序blood<-data.frame(X1=c(76.0,91.5,85.5,82.5,79.0,80.5,74.5,79.0,85.0,76.5,82.0,95.0,92.5),X2=c(50,20,20,30,30,50,60,50,40,55,40,40,20),Y=c(120,141,124,126,117,125,123,125,132,123,132,155,147))#建立數(shù)據框lm.sol<-lm(Y~X1+X2,data=blood)

#進行回歸分析summary(lm.sol)

#匯總分析結果Y=-62.96+2.136X1+0.4002X2.預測：X=(80,40)時，相應Y的概率為0.95的預測區(qū)間.new<-data.frame(X1=c(80,75),X2=c(40,38))lm.pred<-predict(lm.sol,new,interval="prediction",level=0.95)lm.pred求解程序blood<-data.frame(X1=c(76回歸分析的匯總結果Call:lm(formula=Y~X1+X2,data=blood)Residuals:Min1QMedian3QMax-4.0404-1.01830.46400.69084.3274Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)-62.9633616.99976-3.7040.004083**X12.136560.1753412.1852.53e-07***X20.400220.083214.8100.000713***---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1Residualstandarderror:2.854on10degreesoffreedomMultipleR-squared:0.9461,AdjustedR-squared:0.9354F-statistic:87.84on2and10DF,p-value:4.531e-07預測結果如下：

fit

lwr

upr1

123.9699

117.2889

130.6509回歸分析的匯總結果Call:回歸診斷par(mfrow=c(2,2))

#設置畫圖為2x2的格式plot(lm.sol,which=c(1:4))

#模型檢驗4張圖，包括殘差圖、QQ圖和Cook距離圖數(shù)據太少，上面診斷結果并不理想。library(car)#載入程序包Car,vif()函數(shù)在其內round(vif(lm.sol),2)

#計算模型的方差膨脹因子，用2位小數(shù)點的格式展示各變量的方差膨脹因子情況如下：X1X21.961.96可以看到所有參數(shù)估計的VIFj=1/(1-Rj2)值都遠遠小于10，并且接近1。因此這里我們不用擔心多重共線性的問題。回歸診斷par(mfrow=c(2,2))#設置畫圖為2x二項選擇模型當我們考慮多個連續(xù)解釋變量對某個取0-1值的響應變量的影響時，R中常用probit或logit回歸來分析。probit:

-1(P{Y=1})=0+Xlogit:

logit(P{Y=1})=log(P{Y=1}/(1-P{Y=1}))=0+X對二項選擇的probit/logit回歸,R軟件可用glm()處理.fm<-glm(formula,family=binomial(link=probit),data=data.frame)fm<-glm(formula,family=binomial(link=logit),data=data.frame)在用glm()函數(shù)處理二項選擇模型時，公式中響應變量y的輸入形式有兩種：1、y中第一列為對應自變量的響應次數(shù)，第2列是不響應的次數(shù)；2、y是只由0、1構成的向量，分別表示對應自變量取值是不響應還是相應。二項選擇模型當我們考慮多個連續(xù)解釋變量對某個取0-1值的響二項選擇案例1研究小電流對農場動物的影響.選擇了7頭牛,6種電擊強度0,1,3,4,5毫安.給出每種電擊強度70次試驗中牛發(fā)生響應的總次數(shù).試分析電擊對牛的影響。二項選擇案例1研究小電流對農場動物的影響.選擇了7頭牛,6種案例1的程序norell<-data.frame(x=0:5,n=rep(70,6),success=c(0,9,21,47,60,63))norell$Ymat<-cbind(norell$success,norell$n-norell$success)glm.sol<-glm(Ymat~x,family=binomial,data=norell)#logit回歸#glm.sol<-glm(Ymat~x,family=binomial(link=probit),data=norell)summary(glm.sol)預測:pre<-predict(glm.sol,data.frame(x=3.5))p<-exp(pre)/(1+exp(pre));pd<-seq(0,5,len=100)pre<-predict(glm.sol,data.frame(x=d))p<-exp(pre)/(1+exp(pre))norell$y<-norell$success/norell$nplot(norell$x,norell$y);lines(d,p)#points(x,y)、lines(x,y)可以在已有圖中添加點、線案例1的程序norell<-data.frame(x=0:5二項選擇案例250位急性林巴細胞性白血病病人，在入院治療時取得了外轅血中的細胞數(shù)X1、林巴結浸潤等級X2(分為0,1,2,3級);出院后有無鞏固治療X3("1”表示有鞏固治療，"0”表示無鞏固治療).并取得病人的生存時間,Y=0表示生存時間在1年以內，Y=1表示生存時間在1年或1年以上.試分析病人生存時間長短的概率與X1,X2,X3的關系.二項選擇案例250位急性林巴細胞性白血病病人，在入院治療時取案例2的程序life<-data.frame(X1=c(2.5,173,119,10,502,4,14.4,2,40,6.6,21.4,2.8,2.5,6,3.5,62.2,10.8,21.6,2,3.4,5.1,2.4,1.7,1.1,12.8,1.2,3.5,39.7,62.4,2.4,34.7,28.4,0.9,30.6,5.8,6.1,2.7,4.7,128,35,2,8.5,2,2,4.3,244.8,4,5.1,32,1.4),X2=rep(c(0,2,0,2,0,2,0,2,0,2,0,2,0,2,0,2,0,2,0,2,0,2,0),c(1,4,2,2,1,1,8,1,5,1,5,1,1,1,2,1,1,1,3,1,2,1,4)),X3=rep(c(0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1),c(6,1,3,1,3,1,1,5,1,3,7,1,1,3,1,1,2,9)),Y=rep(c(0,1,0,1),c(15,10,15,10)))glm.sol<-glm(Y~X1+X2+X3,family=binomial,data=life)summary(glm.sol)案例2的程序life<-data.frame(多項Logistic回歸模型問題：因變量為多分類變量,自變量為定距/定比變量或虛擬變量。假設因變量有k+1個類別，那么選定一個(如第一個)為參照類，其余各類與參照類相對照，結果可以得到logistic函數(shù)：發(fā)生比：屬于i類與屬于0類的概率之比處理多項Logistic回歸的R指令通常調用格式： library(mlogit) mlogit(formula,data,weights,start=NULL,reflevel=NULL,probit=FALSE,panel=FALSE,shape="long",...)默認數(shù)據格式是”long”,而實際數(shù)據通常是”wide”格式。多項Logistic回歸模型問題：因變量為多分類變量,自變量多分類Logistic回歸示例92年美國總統(tǒng)選舉數(shù)據,內有:投票結果pres92,性別sex,學歷degree,年齡age,年齡段agecat,受教育程度educ.試考察對投票結果有影響的因素。R程序如下：tt=read.csv("e:/tmp/pres92.csv",header=T);head(tt)tt$degree=factor(tt$degree,levels=c("graduatedegree","lthighschool","highschool","juniorcollege","bachelor"))sol=mlogit(pres92~0|age+degree+sex,tt,shape="wide",reflevel="Clinton");summary(sol)sol=mlogit(pres92~0|degree+sex,tt,shape="wide",reflevel="Clinton");summary(sol)#似然比卡方改變較大,age不應刪！多分類Logistic回歸示例92年美國總統(tǒng)選舉數(shù)據,內有:對發(fā)生比率的解讀自變量的發(fā)生比率exp(

)可以這樣表述：定距/比層次的自變量：如果用因變量b類和a類(參照類)相比，自變量每增加一個單位，選擇b類的發(fā)生比是原來的發(fā)生比的exp(

)倍。定類/序層次的自變量：如果因變量b類和a類（參照類）相比，自變量中的某一類選擇b類的發(fā)生比是自變量中的參照類選擇b類的發(fā)生比的exp(

)倍對發(fā)生比率的解讀自變量的發(fā)生比率exp()可以這樣表述：定序回歸模型當我們考察多個連續(xù)解釋變量對某個取有序離散值的響應變量的影響時，適用定序回歸，R中常用MASS包中polr函數(shù)來實現(xiàn)。Probit定序回歸:

-1(P{Yk})=k-

kXLogit定序回歸:

logit(P{Yk})=k-

kX

其中：logit(p)=log[p/(1-p)]polr(formula,data,weights,start,...,subset,na.action,contrasts=NULL,Hess=FALSE,model=TRUE,method=c("logistic","probit","cloglog","cauchit"))定序回歸模型當我們考察多個連續(xù)解釋變量對某個取有序離散值的響定序回歸示例選取商務手機廠商較多考慮的六個功能，然后加上品牌(只涉及諾基亞、摩托羅拉、三星和波導四個品牌)共七個要素，構成我們要研究的影響消費者偏好的要素。rm(list=ls())#清空當前工作空間a=read.csv("D:/phone.csv")#讀入csv格式的訓練數(shù)據attach(a)#將a中的各變量放入工作空間，便于直接調用定序回歸示例選取商務手機廠商較多考慮的六個功能，然后加上品牌示例：1.描述性分析品牌影響、數(shù)碼相機影響…xtabs(~score+W1)#根據score和W1的生成列聯(lián)表xtabs(~score+W2)plot(c(1,5),c(0,1),type=“n”,xlab=“score”,ylab=“Percentage”,main=“W2”)#生成畫圖框架,給出x、y軸標簽標題points(c(1:5),tapply(W2,score,mean),type=“b”) #據score的取值計算的W2的均值:tapply(X,INDEX,FUN=NULL);畫點圖par(mfrow=c(2,2))#生成2x2的圖形plot(c(1,5),c(0,1),type="n",xlab="score",ylab="Percentage",main=“W3")points(tapply(W3,score,mean),type="b")plot(c(1,