用R语言做非参数和半参数回归笔记学习资料

1、用R语言做非参数和半参数回归笔记由詹鹏整理，仅供交流和学习根据南京财经大学统计系孙瑞博副教授的课件修改，在此感谢孙老师的辛勤付出！教材为： Luke Keele: Semiparametric Regression for the Social Sciences. JohnWiley & Sons, Ltd. 2008.第_章 introduction: Global versus Local Statistic一、主要参考书目及说明1、Hardle(1994). Applied Nonparametic Regresstion.较早的经典书2、Hardle etc (2004). Nonpa

2、rametric and semiparametric models: an introduction.Springer.结构清晰3、Li and Racine(2007). Nonparametric econometrics: Theory and Practice.Princeton.较全面和深入的介绍，偏难4、Pagan and Ullah (1999). Nonparametric Econometrics.经典5、Yatchew(2003). Semiparametric Regression for the Applied Econometrician. 例 子不错6、高铁梅(2

3、009).计量经济分析方法与建模：EVIEWS应用及实例(第二版).清华大学出版社.(P127/143)7、李雪松(2008).高级计量经济学.中国社会科学出版社.(P45 ch3)8、 陈强(2010).高级计量经济学及Stata应用.高教出版社.(ch23/24)【其他参看原ppt第一章】二、内容简介方法：移动平均(moving average)核光滑(Kernel smoothing )K近邻光滑(K-NN)局部多项式回归(Local Polynormal)Loesss and Lowess样条光滑(Smoothing Spline )B-splineFriedman Supersmoo

4、ther模型：非参数密度估计非参数回归模型非参数回归模型时间序列的半参数模型Panel data的半参数模型Quantile Regression三、不同的模型形式1、线性模型 linear models2、Nonlinear in variables3、Nonlinear in parametersP小+炊四、数据转换Power transformation （对参数方法）In the GLM framework, models are equally prone （倾向于）to some misspecification （不规范）from an incorrect functional

5、form.It would be prudent （谨慎的）to test that the effect of any independent variable of a model does not have a nonlinear effect. If it does have a nonlinear effect, analysts in the social science usually rely on Power Transformations to address nonlinearity. ADD:检验方法见 Sanford Weisberg. Applied Linear

6、Regression （Third Edition）. A John Wiley & Sons, Inc., Publication.（本科的应用回归分析课教材）ln（为 Jl = O第二章 Nonparametric Density Estimation非参数密度估计直方图Hiatogram1、直方图 Hiatogram2、Kernel density estimate、三种方法2h2h3、K nearest-neighbors estimate二、Histogram对直方图的一个数值解释Suppose x1,xN f（x）, the density function f（x） is unk

7、nown.dx:One can use the following function to estimate f（x）【与x的距离小于h的所有点的个数】三、Kernel density estimateBandwidth: h; Window width: 2h.1、Kernel function 的条件The kernel function K(.) is a continuous function, symmetric(对称的)around zero, that integrates(积分)to unity and satisfies additional bounded conditio

8、ns:K() is symmetric around 0 and is continuous;网泌=1 J颈泌=。诉例feSgr(a) K(z)=0 if |z|=z0 for z0Or(b) IzlK(z) a0 as(4), where is a constant.2、主要函数形式K-tnitlKenwlK(2 16IBiirorni nr box or icctiingukir)| x nui L3510Tnaiitiuki! I osI(1 - I明 乂 Kll 1)Epiin口丽氐it (or qumlrmib4(l f 气 X l(|z| 蒜u 亍时x m：i n2.0362Tri

9、weilii券（i笋项。窿i（ui ij2.3122器（1 一以平X III斜r77MFunh-rttM m F (2 ?r 厂心 m pt 二上 / 2)Fourih*?rikr quariiL 10c +7二x 1( z w 1)3、置信区间f(x)-Bias(x) 1 .96 j 佥了M注其中， 洞六可+； 5对成叽用皿可我祈二:f广（动央血Var如见阪）山。（4、窗宽的选择实际应用中，& = 13收 L3E.W）。其中，s是样本标准差，iqr 是样本分位数级差（interquartile range）四、K nearest-neighbors estimate五、R语言部分da - r

10、ead.table(PSID.txt,header=TRUE)Ihwage - da$lhwage#* bandwidth相等，核函数不同*den1 - density (lhwage,bw=0.45,kernel,epan)den2 - density(lhwage,bw=0.45,kernel=gauss)den3 - density(lhwage,bw=0.45,kernel=biwe)den4 - density(lhwage,bw=0.45,kernel=rect)plot(den4,lty=4,main= ,xlab=Log Hourly Wage”,ylab=Kernel den

11、sity estimates)lines(den3,lty=3,col = red)lines(den2,lty=2, col=green)lines(den1,lty=1,col=blue)#* bandwidth不相等，核函数也不同*den5 - density(lhwage,bw=0.545,kernel=epan)den6 - density(lhwage,bw=0.246,kernel=gauss)den7 - density(lhwage,bw=0.646,kernel=biwe)den8 - density(lhwage,bw=0.214,kernel=rect)plot(den

12、8,lty=4,main= ,xlab=Log Hourly Wage,ylab=Kernel density estimates)lines(den7,lty=3,col = red)lines(den6,lty=2, col=green)lines(den5,lty=1,col=blue)第三章 smoothing and local regression一、简单光滑估计法Simple Smoothing1、Local Averaging 局部均值按照x排序，将样本分成若干部分(intervals or “bins将每部分x对应的y 值的均值作为f(x)的估计。三种不同方法：相同的宽度(e

13、qual width bins) : uniformly distributed.相同的观察值个数(equal no. of observations bins)： k-nearest neighbor.(3)移动平均(moving average)K-NN：等窗宽：2、kernel smoothing 核光滑其中二、局部多项式估计 Local Polynomial Regression1、主要结构局部多项式估计是核光滑的扩展，也是基于局部加权均值构造。local constant regressionlocal linear regressionlowess (Cleveland, 1979

14、)loess (Cleveland, 1988)【本部分可参考:Takezana(2006). Introduction to Nonparametric Regression.(P185 3.7 and P195 3.9)Chambers and Hastie(1993). Statistical models in S. (P312 ch8)】2、方法思路对于每个xi,以该点为中心，按照预定宽度构造一个区间； 在每个结点区域内，采用加权最小二乘法(WLS)估计其参数，并用得 到的模型估计该结点对应的x值对应y值，作为y|xi的估计值(只要这一个点 的估计值)；估计下一个点xj；将每个y|x

15、i的估计值连接起来。【R操作library(KernSmooth) #函数 locpoly()library(locpol) #locpol(); locCteSmootherC()library(locfit) #locfit()#weight funciton: kernel=”tcub”. And“rect”， “trwt”， “tria”， “epan”， “bis(3、每个方法对应的估计形式变量个数 p=0, local constant regression (kernel smoothing)Qlffs)=ix件-禹亍=mina * 每i-i变量个数 p=1, local lin

16、ear regressionminQlffsM莺平罚Cvf 一妍Lowess (Local Weighted scatterplot smoothing) p=i：，=A+A*0A.A)E-wmin口 I /【还有个加权修正的过程，这里略，详见原书或者PPT】Loess (Local regression)p=1,2/=A+A+A+A/9=力4压-片-用u-为成min口 -【还有个加权修正的过程，这里略，详见原书或者PPT】Friedman supersmoothersymmetric k-NN , using local linear fit,varying span, which is d

17、etermined by local CV,not robust to outliers, fast to computesupsmu( ) in R三、模型选择需要选择的内容：(1)窗宽the span；(2)多项式的度the degree of polynomial for the local regression models；(3)权重函数 the weight functions。【其他略】四、R语言部分library(foreign)library(SemiPar)library(mgcv)jacob - read.table(jacob.txt,header=TRUE)#第一部分

18、，简单的光滑估计#1、Kernel Density Estimation#Illustration of Kernel Concepts#Defining the Window Widthattach(jacob)x0 - sort(perotvote)75diffs - abs(perotvote - x0)which.diff - sort(diffs)120#Applying the Tricube Weight#.Tricube functiontricube - function(z) ifelse (abs(z) 1, (1 - (abs(z)A3)A3, 0)#.a - seq(

19、0,1, by=.1)tricube(a)#Figure 2.5plot(range(perotvote), c(0,1), xlab=Perot Vote (%), ylab=Tricube Weight, type=n, bty=l)abline(v=c(x0-which.diff, x0+which.diff), lty=2)abline(v=x0)xwts - seq(x0-which.diff, x0+which.diff, len=250)lines(xwts, tricube(xwts-x0)/which.diff), lty=1, lwd=1)points(x.n, tricu

20、be(x.n - x0)/which.diff), cex=1)#2、Kernel Smoothing#Figure 2.6par(mfrow=c(3,1)plot(perotvote, chal.vote, pch=.”, cex=1.95,xlab=Perot Vote (%), ylab=Challengers Vote Share (%),main=Bandwidth = 4, bty=l)lines(ksmooth(perotvote, chal.vote, bandwidth=4)plot(perotvote, chal.vote, pch=.”, cex=.65,xlab=Per

21、ot Vote (%), ylab=Challengers Vote Share (%),main=Bandwidth = 8”, bty=l)lines(ksmooth(perotvote, chal.vote, kernel=box”,bandwidth=8”), lty=1)plot(perotvote, chal.vote, pch=.”, cex=.65,xlab=Perot Vote (%), ylab=Challengers Vote Share (%),main=Bandwidth = 12”, bty=l)lines(ksmooth(perotvote, chal.vote,

22、 bandwidth=12”), lty=1)* Kernel smoothing中选取box和normal核函数的比较，带宽相等plot(perotvote, chal.vote, pch=.”, cex=.65, xlab=Perot Vote (%), ylab=Challengers Vote Share(%), main=Bandwidth = 8”, bty=l)lines(ksmooth(perotvote, chal.vote, kernel=box,bandwidth=8), lty=1)lines(ksmooth(perotvote, chal.vote, kernel=n

23、ormal,bandwidth=8), lty=2,col=red)#第二部分，LPR模型 #Data Prep For Local Average Regression Step-by-Stepcong - as.data.frame(jacob,2:3)cong - congorder(cong$perotvote),1:2y - as.matrix(cong$chal.vote)x - as.matrix(cong$perotvote)n - length(y)#.tricube - function(z) ifelse (abs(z) 1, (1 - (abs(z)A3)A3, 0)#

24、.x0 - x75diffs - abs(x - x0)which.diff - sort(diffs)120n - xdiffs= which.diffy.n - ydiffs = which.diffweigh=tricube(x.n-x0)/which.diff)mod which.diff, ydiffs which.diff, pch=16, cex=1, col=gray(.80)points(xdiffs = which.diff, ydiffs = which.diff, cex=.85)abline(mod, lwd=2, col=1)text(27.5, 50, expre

25、ssion(paste(Fitted Value of y at , x0) #这里 expression 的用法比较有意思 arrows(25, 47, 15, 37, code =2, length = .10)#2、Now Putting It Together For Local Regression Demonstration.#OLS Fit for Comparisonols - lm(chal.vote perotvote, data=jacob)#The loess fitmodel.loess - loess(chal.vote perotvote, data=jacob,

26、 span = 0.5)#* 默认设置 degree=2, family=gauss, tricube 加权 *n - length(chal.vote)loess - seq(min(perotvote), max(perotvote), length=n)y.loess - predict(model.loess, data.frame(perotvote=x.loess) #得到预测值便于比较#The lowess fitmodel.lowess - lowess(chal.vote perotvote, data=jacob, f = 0.5)#* 默认设置 robust linear

27、 tricube 加权 *n - length(chal.vote)lowess - seq(min(perotvote), max(perotvote), length=n)y.lowess - predict(model.lowess, data.frame(perotvote=x.lowess) #得到预测值便于比较#Figure 2.8plot(perotvote, chal.vote, pch=.,ylab=Challengers Vote Share (%), xlab=Vote for Perot (%), bty=l)lines(x.loess, y.loess)lines(x

28、.lowess, y.lowess)abline(ols)legend(15,20, c(Loess,Lowess, OLS) , lty=c(1,2,1), bty=n, cex=.8)#3、lowess中不同robust的比较m1.lowess - lowess(perotvote, chal.vote, f = 0.5, iter=0)#*没有进行第二步的robust加权估计*m2.lowess - lowess(perotvote, chal.vote, f = 0.5)#*默认iter=3，要进行3次robust加权估计*m0.loess - loess(chal.vote pero

29、tvote, data=jacob,span = 0.5, degree=1, family=symm, iterations=1) #* no robustm1.loess - loess(chal.vote perotvote, data=jacob,span = 0.5, degree=1) #*没有进行第二步的robust加权估计*m2.loess - loess(chal.vote perotvote, data=jacob,span = 0.5, degree=1, family=symm, iterations=3)#*进行3次robust加权估计*plot(perotvote,

30、 chal.vote, pch=., ylab=Challengers Vote Share (%), xlab=Vote for Perot (%)lines(m1.lowess)lines(sort(perotvote), m1.loess$fitorder(perotvote), lty=3, col=green)lines(sort(perotvote), m0.loess$fitorder(perotvote), lty=9,col=18)lines(m2.lowess, lty=2, col=red)lines(sort(perotvote), m2.loess$fitorder(

31、perotvote), lty=4, col=blue)第四章样条估计spline一、基本思想按照x将样本分成多个区间，对每个区间分别进行估计。不同于核估计，这里 不用移动计算，从而减小了计算量。二、最简单的形式Linear Spline with k knots:二%+ tO+A； CO+，AAhCO+其中，”弘r o 皿，)=cx-cx三、其他样条模型1、p次样条估计二次样条 Quadratic Spline (basis functions with k knots)+ ope* 虹HM (x-rR 心三次样条 Cubic Spline (with k knots, use quadra

32、tic basis functions)p = % + 4工，虹 * * +-就+，用工-r) j * Ep-order spline (with k knots)p = % + W 4用A 花+4威Ct-R片槌2、B-splines (with k knots cubic B-spline basis)珀沪十牌3=%jc? +其中，*)= 2*祯,5=30 othEiwise畦3 =二A尸与 5 * AgCnCjCnl3、Natural Splines以上估计方法对结点(knots)之间的估计比较准确，但对边界的拟合效果较 差。自然样条的思想是，在自变量最小值和最大值处各增加一个结点，用线性

33、 模型拟合边界位置的样本点。4、k的选择和模型比较采用AIC准则四、光滑样条 smoothing spline基于如果目标得到参数估计值min口五、模型比较的F检验IF _ g。-河使 5六、R语言部分library(foreign)jacob - read.dta(jacob.dta)attach(jacob)#第一部分，B样条和natural B样条library(splines)#* P61 Perform Spline Regression *m.bsp - lm(chal.votebs(perotvote, df=5), data=jacob)#* 3次B样条公式:df=k+3 (不

34、含常数项)m.nsp - lm(chal.votens(perotvote, df=5), data=jacob)#* df=5对应结点为4个；3次natural B样条公式：df=k+1perot - seq(min(perotvote), max(perotvote), length=312)bsfit - predict(m.bsp, data.frame(perotvote=perot)nsfit - predict(m.nsp, data.frame(perotvote=perot)AlC(m.bsp) #计算 AIC 值#第二部分，光滑样条估计plot(perotvote, cha

35、l.vote, pch=.”,ylab=Challengers Vote Share (%), xlab=Vote for Perot (%),bty=l”, main = df = 2, cex.main = .95)lines(smoothspline(perotvote, chal.vote, df=2)#厢#第三部分，置信区间library(splines)m.nsp - lm(chal.votens(perotvote, df=4), data=jacob)perot - seq(min(perotvote), max(perotvote), length=312)nsfit - p

36、redict(m.nsp, inteval=confidence, se.fit=TRUE,data.frame(perotvote=perot)#Figure 3.8plot(perotvote, chal.vote, pch=., ylab=Challengers Vote Share (%), xlab=Vote for Perot (%)lines(perot, nsfit$fit)lines(perot, nsfit$fit + 1.96*nsfit$se.fit, lty=2)lines(perot, nsfit$fit - 1.96*nsfit$se.fit, lty=2)#*偏

37、移调整的置信区间#* Figure 3.9#* Overlay Natural cubic B-pline Fit and Confidence Intervalslibrary(mgcv)m.smsp - gam(chal.vote s(perotvote, bs=cr, k=4, fx=TRUE)plot(m.smsp, rug=FALSE, se=TRUE,ylab=Challengers Vote Share (%), xlab=Vote for Perot (%),residual=FALSE, shift=33.88) #绘制调整以后的上下界lines(perot, nsfit$f

38、it + 1.96*nsfit$se.fit, lty=3)lines(perot, nsfit$fit - 1.96*nsfit$se.fit, lty=3)legend(3,47, c(natural cubic B-spline Pointwise Bands, Bias Adjusted Bands), lty=c(3, 2), bty=n)#第四部分，模型比较ols - lm(chal.vote 1)#Use Automated R F-test Function Anovaanova(ols, m.nsp)anova(ols, m.smsp)#*三个不同的函数做smoothing

39、spline，进行比较plot(x,y, xlab=X”, ylab=Y)#方法1library(SemiPar)fit - spm(y f(x)lines(fit, se=FALSE, lwd=1,lty=1,col=1)#方法2library(pspline)lines(sm.spline(x,y, df=31),lty=5,col=5)#方法3lines(smooth,spline(x,y,df=31),lty=6,col=6)legend(0,-0.8,c(spm”,sm.spline”,smooth.spline),lty=c(1,5,6),col=c(1,5,6),cex=0.8,

40、bty=n)第五章 Automated Smoothing Techniques 自动光滑技术一、Span by Cross-ValidationIGCV*) =瞄&)=也-孑5(n-40口其中s指窗宽span。【CV和GCV在LPR中表现不佳】二、自动光滑技术两种方法：1、采用MLE (极大似然估计)；2、采用CV选择。1、MLE方法/(V) = -inhi(2j+ln|V| +Yy得到。其中，是随机效应(the random effect)的方差估计，氏是随机误差项(the error term)的方差估计2、最小化CV或GCVGCV(A)=n-n-v(sjfn-n-hrtSjf三、R语言

41、部分setwd(D:/ 课程 /nonparameter regression/2010/ch3”)jacob - read.table(jacob.txt,header=TRUE)attach(jacob)#第一部分，计算GCV并绘图library(locfit)alpha - seq(0.2,0.8, by=0.05)plot(gcvplot(chal.voteperotvote, data=jacob, alpha=alpha), type=o)#第二部分，比较MLE方法和robust得到的估计结果* P89 Figure 4.2 *plot(perotvote, chal.vote,

42、pch=.”, cex=1.75,ylab=Challengers Vote Share (%), xlab=Vote for Perot (%)fit - locfit(chal.voteperotvote, data=jacob, alpha=0.2)lines(fit)lines(lowess(perotvote, chal.vote, f = 0.2), lty=2,col=2) # robustlines(lowess(perotvote, chal.vote, f = 0.2,iter=0), lty=3,col=3) # no robustlegend(23,16, c(locf

43、it with likelihood”,robust lowess”,norobust lowess), lty=1:3,col=1:3,bty=n”,cex=0.8)#第三部分，在光滑样条中是哟花姑娘GCVlibrary(pspline) #使用其函数 sm.spline()plot(perotvote, chal.vote, pch=.”,ylab=Challengers Vote Share (%), xlab=Vote for Perot (%),main = cubic smoothing spline With GCV)lines(sm.spline(perotvote, chal

44、.vote, cv=F),lty=2,col=2) #GCV Smoothing Selectionlines(smooth,spline(perotvote, chal.vote, cv=F),lty=3,col=3)#第四部分，基于GCV的LS方法和MLE方法比较library(mgcv) #使用其函数 gam()library(SemiPar) #使用其函数 spm()smsp1 - gam(chal.vote s(perotvote, bs=cr) # 用 LS 方法估计，GCVsmsp2 - spm(chal.vote f(perotvote) #光滑样条的混合模型表示，用极大似然方

45、法估计！smsp1 # or summary(smsp1)summary(smsp2)* P92 Figure 4.3 *par(mfrow=c(1,2)plot(smsp1, rug=FALSE, se=FALSE,ylab=Challengers Vote Share (%), xlab=Vote for Perot (%),residual=TRUE, shift=33.88, bty=l, main=GCV Smoothing)points(perotvote, chal.vote, pch=., cex=1.75)plot(perotvote, chal.vote, pch=., c

46、ex=1.75,ylab=Challengers Vote Share (%), xlab= Vote for Perot (%), main=Likelihood Smoothing, bty=l)lines(smsp2, se=FALSE, lwd=1, rug=FALSE)第六章 Additive and Semiparametric Regression Models可加回归模型和半参数回归模型十KChH 十尢伏）十用褊十十月*R语言部分 jacob - read.table(jacob.txt, header=TRUE)attach(jacob)#第一部分，广义可加模型library

47、(mgcv)gam1 - gam(chal.vote s(perotvote, bs=cr) + s(checks, bs=cr), data=jacob) summary(gam1)#画图par(mfrow = c(1,2)plot(gam1, select=1, rug=F, se=TRUE,ylab=Challengers Vote Share (%), xlab=Vote for Perot (%), bty=l”, shift=36.6672)points(perotvote, chal.vote, pch=., cex=1.75)plot(gam1, select=2, rug=F

48、, se=TRUE,ylab=Challengers Vote Share (%), xlab=Number of Overdrafts, bty=l”, shift=36.6672)points(checks, chal.vote, pch=., cex=1.75)#进行检查gam.check(gam1)#模型比较的卡方检验#OLS Modelols1 - gam(chal.vote perotvote + checks, data=jacob)#或 ols - lm(chal.vote perotvote + checks, data=jacob) # 速度更快 summary(gam1)

49、$n #自由度或者用 sum(gam1$edf)#Chi sqaured test：将线性模型与非参的可加模型进行比较 #deviance计算模型的变异度LR - summary(gam1)$n*(log(deviance(ols1) - log(deviance(gam1)df - sum(gam1$edf) - sum(ols1$edf)1 - pchisq(LR, df)#第二种比较方法，采用anovaanova(ols1, gam1, test=Chisq)#第二部分，半参数模型library(mgcv)#* P123 Baseline Model * no transformatio

50、ns *ols - gam(chal.vote exp.chal + chal.spend + inc.spend + pres.vote + checks +marginal + partisan.redist + perotvote, data=jacob)* Test each continuous covariategam1 - gam(chal.vote exp.chal + s(chal.spend, bs=cr) + inc.spend + pres.vote+ checks + marginal + partisan.redist + perotvote, data=jacob

51、)gam2 - gam(chal.vote exp.chal + chal.spend + s(inc.spend, bs=cr) + pres.vote+ checks + marginal + partisan.redist + perotvote, data=jacob)gam3 - gam(chal.vote exp.chal + chal.spend + inc.spend + s(pres.vote, bs=cr)+ checks + marginal + partisan.redist + perotvote, data=jacob)gam4 - gam(chal.vote ex

52、p.chal + chal.spend + inc.spend + pres.vote + s(checks, bs=cr)+ marginal + partisan.redist + perotvote, data=jacob)gam5 - gam(chal.vote exp.chal + chal.spend + inc.spend + pres.vote + checks + marginal+ partisan.redist + s(perotvote, bs=cr), data=jacob)anova(ols,gam1,test=Chi)anova(ols,gam2,test=Chi

53、)anova(ols,gam3,test=Chi)anova(ols,gam4,test=Chi)anova(ols,gam5,test=Chi)第七章 Generalized Additive Models广义可加模型 一、广义线性模型GLMprobit model: E 舞已尽)；攻传)-垂期Dlogit model: PL 五Q司；漏 1+ #f内0|对=翼疑=1|对=?1possion regression:此 诚二、广义可加模型X.=Xa=K兽=E(F| 与，与 a 一孕)%=_fi号9 一, 力或瑞=e阪十U(%)十-十4X3三、估计方法MLE: use Newton-Raphso

54、n algorithmIRLS: backfitting algorithm (in ch5)四、假设检验LR=-2 (LogLikelihood0 logLikelihoodl)【这是两个模型进行比较】五、R语言部分setwd()library(foreign)war - read.dta(PHdata.dta)attach(war)#第一部分，用广义可加模型做Logit回归【与glm的用法区别】library(mgcv)war.glm - gam(dispute nudem + nugrow + allies + contig + nucapab + trade+ sumdisp + s(

55、year, bs=cr), data = war, family=binomial)summary(war.glm)#进行预测war.smooth7 - gam(dispute s(nudem, bs=cr) + nugrow + allies + contig + nucapab + trade+ sumdisp + s(year, bs=cr), data = war, family=binomial)summary(war.smooth7)new - data.frame(nudem=nudem, nugrow=nugrow, allies=allies, contig=contig,

56、nucapab=nucapab,trade=trade, sumdisp=sumdisp, year=year)sa - predict(war.smooth7, newdata=new) #* eta(i) = XB+f(X)ui - plogis(sa) #* plogis为logistic分布的分布函数，计算概率summary(ui)summary(war.smooth7$fit)sa2 - predict(war.smooth7, newdata=new,type=terms)#*可以得到模型的每一项(terms)的预测值sa2,7 #* f(nudem)那一列#第二部分，用广义可加模

57、型做多分类logit模型(顺序)library(foreign)couples - read.dta(D:/课程/nonparameter regression/2010/Rprog/Chapter6/couples.dta)attach(couples)library(VGAM)gam1 - vgam(violence chabting + minority + s(fagunion) + s(misolatn) + s(ecndisad) + alcdrug+ s(duryrs) + s(econ_alc) + s(disagmnt) + s(comstyle),cumulative(pa

58、rallel=T), data=couples)summary(gam1)par(mfrow=c(3,3)plot(gam1,which.term=3, se=TRUE, rug=FALSE, bty=l)plot(gam1,which.term=4, se=TRUE, rug=FALSE, bty=l)plot(gam1,which.term=5, se=TRUE, rug=FALSE, bty=l)plot(gam1,which.term=7, se=TRUE, rug=FALSE, bty=l)plot(gam1,which.term=8, se=TRUE, rug=FALSE, bty

59、=l)plot(gam1,which.term=9, se=TRUE, rug=FALSE, bty=l)plot(gam1,which.term=10, se=TRUE, rug=FALSE, bty=l)【关于预测方法的说明没有看明白】#第三部分，用广义可加模型做possion模型library(foreign)Scourt - read.dta(scourt.dta)attach(Scourt)library(mgcv)ml - gam(nulls tenure + congress + unified, data=Scourt,family=poisson) summary(ml)m2

60、 - gam(nulls tenure + s(congress, bs=cr) + unified,data=Scourt,family=poisson)summary(m2)anova(m1, m2, test=Chisq)#进行预测predict.data - data.frame(expand.grid(list(congress=seq(1, 101, by=1), tenure =10.35, unified = 0)#Created Fitted Values Using Fake Data 方法 1predict.fit - predict(m2, newdata = pred