多元统计分析习题操作及分析

例2.1>data2.1<-read.table("clipboard",header=T)>lm.salary<-lm(y~x1+x2+x3+x4,data=data2.1)>summary(lm.salary)Call:lm(formula=y~x1+x2+x3+x4,data=data2.1)Residuals:Min1QMedian3QMax-12924.2-4588.1-269.61756.225215.7Coefficients:EstimateStd.ErrortvaluePr(>t|)(Intercept)48386.062011237.28824.3060.000155***x11.68310.130212.9295.01e-14***x2-34.5520130.2602-0.2650.792570x3-13.000413.7882-0.9430.353043x4808.3223547.80171.4760.150144Signif.codes:0‘***'0.001‘**'0.01‘*'0.050.1

Residualstandarderror:7858on31degreesoffreedomMultipleR-squared:0.919,AdjustedR-squared:0.9086F-statistic:87.95on4and31DF,p-value:<2.2e-16从以上输出结果可以看出，回归方程的F值为87.95,相应的P值为，说明回归方程是显著的，但t检验对应的p值则显示常数项和xl是显著的，而x2,x3,x4不显著。变量选择：采用R软件中的step()过程可以完成逐步回归的过程。程序如下：用“一切子集回归法”来进行逐步回归>lm.step<-step(lm.salary,direction="both")Start:AIC=650.4ly~xl+x2+x3+x4-x2l4.3448e+06-x2l4.3448e+06-x3l5.4896e+07<none>-x4ll.3445e+08-xlll.0323e+l0DfSumofSqRSSAICl.9l86e+09648.49l.9692e+09649.43l.9l43e+09650.4l2.0487e+09650.85l.2237e+l07l5.l9Step:AIC=648.49y~x1+x3+x4-x31-x316.2078e+07<none>-x411.3011e+08+x214.3448e+06-x111.0341e+10DfSumofSqRSSAIC1.9807e+09647.641.9186e+09648.492.0487e+09648.851.9143e+09650.411.2259e+10713.26RSSAIC1.9807e+09647.64RSSAIC1.9807e+09647.641.9186e+09648.491.9692e+09649.432.2771e+09650.661.3635e+10715.09+x316.2078e+07+x211.1527e+07-x412.9640e+08-x111.1654e+10<none>Step:AIC=647.64y~x1+x4DfSumofSq采用全部自变量作回归时，AIC=650.41,如果去掉变量x2,AIC值减小为648.49；如果去掉变量x3，AIC值减小为649.43；如果去掉变量x4，AIC值增大为650.85；如果去掉变量x1，AIC值增大为715.09.由于去掉x2,AIC值达到最小，所以R软件去掉x2进入第二轮计算，

此时AIC=648.49，如果去掉变量x3,AIC值减小为647.64；如果去掉其他变量或增加变量，都会使AIC值增大，因此R软件去掉x3进入第三轮计算，此时AIC=647.64，无论去掉那个变量或者增加哪个变量，AIC值都会增大，所以停止计算，得到最优回归模型，即y关于x1和x4的线性回归模型。现在用函数summary(lm.step)来得到回归模型的如下汇总信息：>summary(lm.step)Call:lm(formula=y~x1+x4,data=data2.1)Residuals:Min1QMedian3QMax-13632-4759-615176125076Coefficients:ErrortvaluePr(>|t|)5265.218ErrortvaluePr(>|t|)5265.2187.9953.18e-09***(Intercept)42097.1650.11713.9342.22e-15***467.6712.2220.0332*x11.631x41039.260Signif.codes:0‘***'Signif.codes:0‘***'0.001‘**'0.01‘*'0.050.1‘'1Residualstandarderror:7747on33degreesoffreedomMultipleR-squared:0.9162,AdjustedR-squared:0.9111F-statistic:180.4on2and33DF,p-value:<2.2e-16注意到常数项、x1和x4都是显著的，模型也是显著的，所以可以得到最优回归方程Y=回归诊断：分别采用residuals(),rstandard()和rstudent()来计算普通残差、标准化残差和学生化残差，程序如下：y.res<-residuals(lm.salary)y.rst<-rstandard(lm.step)print(y.rst)1234567891011121314-0.174811711.03650457-1.641440643.43509088-0.693888980.21730074-0.28221956-0.57391074-1.138600820.14942833-0.08074416-0.70313518-0.08074416-1.9087926615161718192021222324

25262728252627280.86519102-0.04532806-0.63719129-0.695292950.27405524-0.464535240.16499595-0.071895450.222034020.418839660.25844627-0.12669958-0.01045782-0.294609572930313233343536-0.15151718-0.896744310.474241840.64162319-0.848800140.184959113.32423970-0.62075537>y.fit<-predict(lm.step)>plot(y.res~y.fit)从标准化残差可以看出，4号点和35号点的标准化残差大于3，因此我们认定第4号和第35号观测值是异常点。33从残差图可以看出，残差的分布有随预测值增大而增大的趋势，所以同方差的基本假定可能不成立。尝试采用对数变换来解决方差非齐问题。程序如下：>y.fit<-predict(lm.step)>plot(y.res~y.fit)>plot(y.rst~y.fit)>lm.step_new<-update(lm.step,log(.)~.)>y.rst<-rstandard(lm.step_new)>y.fit<-predict(lm.step_new)>plot(y.rst~y.fit)

.yOO00OO°°□O.yOO00OO°°□Oo00°°O口％oo0DO护O00O0O32101211.211.411.611.812.012.2y.fit比较可以发现，对模型进行对数化变换后残差散点图有所改善，只有35号点是异常点，这里做一个简单的处理，去掉4号和35号观测值，重复上述回归分析和残差分析过程，可以得到新的标准化残差图，发现，残差的分布有了很大的改进，全部落在【-2,2】的带状区域内，上述分析过程的R程序如下：回归诊断：一般的方法ResidualsvsFittedResidualsvsFittedNormalQ-QResidualsvsFittedResidualsvsFittedNormalQ-Q10001o・Fittedvalues54■211■611■811Scale-Location0100■2JIJI4JIJIGJIJI0JIJI10001o・Fittedvalues54■211■611■811Scale-Location0100■2JIJI4JIJIGJIJI0JIJI-2-1012432o1-2-3-TheoreticalQuantilesResidualsvsLeverage432101-2-3-Cook'sdistance36a5a4a3a2a10.50.51FittedvaluesLeverage>par(mfrow=c(2,2))>plot(lm.step_new)>influence.measures(lm.step_new)Influencemeasuresoflm(formula=log(y)~x1+x4,data=data2.1):dfb.1_dfb.x1dfb.x4dffitcov.rcook.dhatinf0.0273880.05593-0.057229-0.084111.1612.42e-030.03730.03730.064620.3337150.07828-0.2770490.352021.0884.11e-020.095330.348810-3.494260.912364-3.860831.2553.93e+000.6054*4-0.5801420.539470.3906961.149420.6493.63e-010.1404*5-0.0705730.033140.018032-0.150681.0647.64e-030.03246-0.034050-0.046230.0566280.073131.1831.83e-030.078770.0006310.00052-0.000879-0.001081.1954.01e-070.08238-0.028534-0.008120.017639-0.050961.1248.91e-040.031790.1011010.18171-0.205807-0.313530.9923.21e-020.0542100.022721-0.00740-0.0073800.046531.1267.43e-040.031711-0.0111970.006970.002037-0.025041.1332.15e-040.0335120.1389610.20784-0.237048-0.303321.1043.07e-020.0891-0.0111970.006970.002037-0.025041.1332.15e-040.03350.285691-0.11607-0.238545-0.439931.1676.43e-020.1502150.3624740.13639-0.3227440.385191.0874.90e-020.1036160.0334440.01349-0.0301640.035671.2254.37e-040.105717-0.0882630.12602-0.018114-0.253260.9972.11e-020.0410180.0823540.16817-0.172082-0.252911.0762.14e-020.064619-0.018752-0.002710.014852-0.019881.2081.36e-040.0923200.0102250.00752-0.015035-0.021781.1581.63e-040.053621-0.0033280.00381-0.000231-0.008761.1402.64e-050.038122-0.0255880.02930-0.001775-0.067331.1281.55e-030.0381230.041709-0.00255-0.0188490.079691.1112.17e-030.030924-0.088167-0.062340.1288780.188261.0921.19e-020.0531250.007481-0.008000.0002460.019231.1381.27e-04260.1267650.21456-0.1641900.260581.1432.29e-020.095827-0.133262-0.019280.105546-0.141301.1866.82e-030.0923280.0100030.01053-0.015693-0.021101.1701.53e-040.062929-0.0253570.021650.001794-0.060921.1261.27e-030.035330-0.211050-0.114080.202643-0.231361.2011.82e-020.1197310.2367640.09713-0.2142270.252751.1672.16e-020.106332-0.054779-0.107670.1135420.168521.1219.63e-030.060933-0.0924630.057580.016818-0.206761.0131.41e-020.033534-0.0857860.134220.0425390.227721.2181.76e-020.128735-0.437405-0.282590.6311170.939190.3412.02e-010.0520*36-0.0531400.010070.020734-0.104861.0963.74e-030.0312上图给出了逐步回归模型lm.step-new的四个回归诊断图，从这四个途中可以看出，除了第3,4，35号观测值，残差-拟合图中的点基本1151115111511151上呈随机分布模式；正态Q-Q图中的点基本落在直线上，表明残差服从正态分布；大小-位置图和残差-杠杆图以小组的形式存在并且离中心不远，这说明第3，4,35号观测值可能是异常点和强影响点。influence.measures(lm.step_new)给出了几个诊断统计量，注意到第3，4,35号观测值的右端有一个星号，说明第3，4,35号观测值被诊断为影响点。回归预测：preds<-data.frame(x1=20000,x4=20)>predict(lm.step,newdata=preds,interval="prediction",level=0.95)fitlwrupr95493.0978187.28112798.9广义线性模型Logistic模型首先打开数据，然后将数据选中导入软件中建立opinion关于age和sex的logistic回归模型模型汇总，给出模型回归系数的估计和显著性检验等data3.1<-read.table("clipboard",header=T)data3.1xy20110012180301601612213617024160110181251120100151702217016118121170906020116112015190801002212419010018130160130231100>glm.logit<-glm(y~x,family=binomial,data=data3.1)Warningmessage:glm.fit:拟合機率算出来是数值零或一>summary(glm.logit)Call:glm(formula=y~x,family=binomial,data=data3.1)DevianceResiduals:Min1QMedian3QMax-1.21054-0.054980.000000.004331.87356Coefficients:EstimateStd.ErrorzvaluePr(>|z|)(Intercept)-21.280210.5203-2.0230.0431*x1.64290.83311.9720.0486*Signif.codes:0‘***'0.001‘**'0.01‘*'0.050.1‘'1(Dispersionparameterforbinomialfamilytakentobe1)Nulldeviance:62.3610on44degreesoffreedomResidualdeviance:6.1486on43degreesoffreedomAIC:10.149NumberofFisherScoringiterations:9可以看出，回归模型的系数在0.1%的水平上显著，性别和年龄两个指标分别都是显著的，说明不同性别对服务产品的观点显著不同，并且年龄同样如此。回归模型为一=3.49978-1.77815sex-0.07156age对数线性模型采用R软件中的广义线性模型过程glm()来建立泊松对数线性模型首先打开数据，然后将数据选中导入软件中建立y关于x1和x2的泊松对数线性回归模型模型汇总，给出模型回归系数的估计和显著性检验等>data3.2<-read.table("clipboard",header=T)>data3.2Nox1x2x3y1111310142211300143362501144836013556622055

66272902277123101288524209599233702210101028033111152360661212332403013131823016141442360421515872605916165026016171718280618181113101231919183201520202021016212112290142222921014232317320132424282503025255530014326269400627271019010282829293030313132323333343435353636373738383939404041414242434344444545464647474848494947220537618142383212819201710301131918119242411131301741435120112711067201244122129728142223161340112463316536211263835139725173626132112513151221302505022321135151412512652523235110535356211705454244111355551632115565622261515757252116585813361059591237110>glm.ln<-glm(y~x1+x2+x3,family=poisson(link=log),data=data3.2)>summary(glm.ln)Call:glm(formula=y~x1+x2+x3,family=poisson(link=log),data=data3.2)DevianceResiduals:Min1QMedian3QMax-6.0569-2.0433-0.93970.792911.0061Coefficients:EstimateStd.ErrorzvaluePr(>z|)(Intercept)1.94882590.135619114.370<2e-16***x10.02265170.000509344.476<2e-16***x20.02274010.00402405.6511.59e-08***x3-0.15270090.0478051-3.1940.0014**Signif.codes:0‘***'0.001‘**'0.01‘*'0.05‘.'0.1‘'1(Dispersionparameterforpoissonfamilytakentobe1)Nulldeviance:2122.73on58degreesoffreedomResidualdeviance:559.44on55degreesoffreedomAIC:850.71NumberofFisherScoringiterations:5可以得到回归模型：In=0.54798+0.02772x1+0.65596x2.从检验结果可以看出家庭年收入(x1)和是否有私家车(x2)的系数都是显著的，说明家庭年收入(x1)和是否有私家车(x2)对家庭一年外出旅游次数(y)有重要影响。家庭年收入的回归系数为0,02772，表明保持其他预测变量不变，收入每增加1万元，一年外出旅游次数对数均值将增加0.2772，是否有私家车的回归系数为0.65596，表明保持其他预测变量不变，有私家车一年外出旅游次数对数均值将增加0.65596.在因变量的初始尺度上解释回归系数，指数化系数为：>exp(coef(glm.ln))(Intercept)x1x2x37.02044031.02291021.02300070.8583864可以看出，表明保持其他预测变量不变，收入每增加1万元，一年外出旅游次数将乘以1.028110，表明保持其他预测变量不变，有私家车一年外出旅游次数将乘以1.926978.聚类分析系统聚类法首先采用系统聚类的最小距离法进行聚类，样品之间的距离采用欧氏距离来度量，将10种葡萄酒看成10类，分别计算各类之间的距离，容易求得6和10之间的距离最小，因此把它们合并为一个新类，记为G,然后采用欧氏距离计算各类之间的距离，发现3和4之间的距离最小，于是把它们合并为一个新类，记为G2;如此一直下去，知道把所有的10种酒合并为一类，聚类合并的顺序如下以上聚类过程的R程序为：数据读入到data4.1中，采用欧氏距离计算相似矩阵d,采用最小距离法(single)聚类，绘制聚类树状图。〉data4.1<-read.table("clipboard",header=T)

>data4.1x1x2x3x4x5x614.654.225.014.504.154.1226.326.116.216.856.526.3334.874.604.954.154.024.1144.884.684.434.124.034.1456.736.656.726.136.516.3667.457.567.607.807.207.1878.108.238.017.958.318.2688.428.548.127.888.267.9896.456.816.526.316.276.06107.507.327.427.527.106.95>d<-dist(data4.1,method="euclidean",diag=T,upper=F,p=2)>HC<-hclust(d,method="single")>plot(HC)ClusterDendrogramthgie34590thgie345901dhclust(*,"single")从图中可以看出，如果取合并距离为4，则10种酒可以分为两类，第一类为6,10,7,8,2,5,9；第二类为1,3,4；如果取合并距离为2则10种酒可以分为三类，第一类为6,10,7,8；第二类为2,5,9；第三类为1,3,4；其次采用系统聚类的最大距离法进行聚类，样品之间的距离采用欧氏距离来度量，聚类的结果和最小距离法的结果一样，只是合并距离有所不同，其程序和结果如下：>HC<-hclust(d,method="complete")>plot(HC)

thgieClusterDendrogram48620thgieClusterDendrogram48620d

hclust(*,"complete")K均值聚类法将eg4.2.xls数据读入到data4.2中，聚类的个数为4,随机集合#的个数为20,算法为"Hartigan-Wong",对分类结果进行排序并查看分类情况data4.2<-read.table("clipboard",header=T)data4.2x1x2x3x4x5x6x7x8北京6905.512265.881923.711562.553521.203306.821523.32975.37天津6663.311754.981763.441174.622699.532116.011415.39836.82河北3927.261425.991372.25809.851526.601203.99955.95387.40山西3558.041461.901327.78832.741487.661419.43851.30415.44内蒙古4962.402514.091418.601162.872003.541812.071239.36765.13辽宁5254.961854.631385.62929.371899.061614.521208.30643.15吉林4252.851769.471468.29839.311541.371468.341108.51562.48黑龙江4348.451681.881185.96723.581363.621190.871082.96476.89上海8905.952053.812225.681826.223808.413746.381140.821394.86江苏6060.911772.061187.741193.812262.192695.52962.45647.06浙江7066.222138.991518.061109.423728.232816.121248.90811.51安徽5246.761371.011501.39690.661365.011631.28

907.58467.77福建6534.941494.961661.841179.8