用R语言做非参数和半参数回归笔记

由詹鹏整理，仅供交流和学习

根据南京财经大学统计系孙瑞博副教授的课件修改，在此感谢孙老师的辛勤付出！ 教材为：Luke Keele: Semiparametric Regression for the Social Sciences. John Wiley & Sons, Ltd. 2008.

------------------------------------------------------------------------- 第一章 introduction: Global versus Local Statistic 一、主要参考书目及说明

1、Hardle(1994). Applied Nonparametic Regresstion. 较早的经典书

2、Hardle etc (2004). Nonparametric 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、高铁梅（2009）. 计量经济分析方法与建模：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-spline

——Friedman Supersmoother 模型：

——非参数密度估计 ——非参数回归模型 ——非参数回归模型

——时间序列的半参数模型 ——Panel data 的半参数模型 ——Quantile Regression 三、不同的模型形式 1、线性模型linear models

2、Nonlinear in variables

3、Nonlinear in parameters

四、数据转换 Power transformation（对参数方法）

In the GLM framework, models are equally prone(倾向于) to some misspecification（不规范） from an incorrect functional 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 Regression (Third Edition). A John Wiley & Sons, Inc., Publication.（本科的应用回归分析课教材）]

---------------------------------------------------------------------------- 第二章 Nonparametric Density Estimation 非参数密度估计 一、三种方法

1、直方图 Hiatogram 2、Kernel density estimate

3、K nearest-neighbors estimate

二、Histogram 对直方图的一个数值解释

Suppose x1,…xN – f(x), the density function f(x) is unknown.

One can use the following function to estimate f(x)

【与x的距离小于h的所有点的个数】 三、Kernel density estimate

Bandwidth: 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 conditions: (1) K() is symmetric around 0 and is continuous; (2)

,

,

;

(3) Either

(a) K(z)=0 if |z|>=z0 for z0 Or

(b) |z|K(z) à0 as (4)

2、主要函数形式

; , where

is a constant.

3、置信区间

其中，

4、窗宽的选择

实际应用中，

位数级差（interquartile range） 四、K nearest-neighbors estimate

。其中，s是样本标准差，iqr是样本分

五、R语言部分

da <- read.table(\lhwage <- da$lhwage #*** bandwidth 相等，核函数不同 *** den1 <- density(lhwage,bw=0.45,kernel=\den2 <- density(lhwage,bw=0.45,kernel=\den3 <- density(lhwage,bw=0.45,kernel=\den4 <- density(lhwage,bw=0.45,kernel=\plot(den4,lty=4,main=\lines(den3,lty=3,col = \lines(den2,lty=2, col=\lines(den1,lty=1,col=\#*** bandwidth 不相等，核函数也不同 *** den5 <- density(lhwage,bw=0.545,kernel=\den6 <- density(lhwage,bw=0.246,kernel=\den7 <- density(lhwage,bw=0.646,kernel=\den8 <- density(lhwage,bw=0.214,kernel=\plot(den8,lty=4,main=\lines(den7,lty=3,col = \lines(den6,lty=2, col=\lines(den5,lty=1,col=\---------------------------------------------------------------------------- 第三章 smoothing and local regression 一、简单光滑估计法 Simple Smoothing 1、Local Averaging 局部均值

按照x排序，将样本分成若干部分（intervals or ―bins‖）；将每部分x对应的y值的均值作为f(x)的估计。 三种不同方法：

（1）相同的宽度（equal width bins）：uniformly distributed.

（2）相同的观察值个数（equal no. of observations bins）：k-nearest neighbor. （3）移动平均（moving average） K-NN：

等窗宽：

2、kernel smoothing 核光滑

其中，

二、局部多项式估计Local Polynomial Regression 1、主要结构

局部多项式估计是核光滑的扩展，也是基于局部加权均值构造。 ——local constant regression ——local linear regression ——lowess (Cleveland, 1979) ——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、方法思路

（1）对于每个xi，以该点为中心，按照预定宽度构造一个区间；

（2）在每个结点区域内，采用加权最小二乘法（WLS）估计其参数，并用得到的模型估计该结点对应的x值对应y值，作为y|xi的估计值（只要这一个点的估计值）； （3）估计下一个点xj；

（4）将每个y|xi的估计值连接起来。 【R操作

library(KernSmooth) #函数locpoly()

library(locpol) #locpol(); locCteSmootherC() library(locfit) #locfit()

#weight funciton: kernel=‖tcub‖. And ―rect‖, ―trwt‖, ―tria‖, ―epan‖, ―bisq‖, ―gauss‖ 】

3、每个方法对应的估计形式

（1）变量个数p=0, local constant regression (kernel smoothing)

min

（2）变量个数p=1, local linear regression

min

（3）Lowess (Local Weighted scatterplot smoothing) p=1:

min

【还有个加权修正的过程，这里略，详见原书或者PPT】 （4）Loess (Local regression) p=1,2:

min

【还有个加权修正的过程，这里略，详见原书或者PPT】 （5）Friedman supersmoother

symmetric k-NN, using local linear fit,

varying span, which is determined by local CV, not robust to outliers, fast to compute supsmu( ) 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(\############################################################################### #第一部分，简单的光滑估计 #1、Kernel Density Estimation #Illustration of Kernel Concepts #Defining the Window Width attach(jacob) x0 <- sort(perotvote)[75] diffs <- abs(perotvote - x0) which.diff <- sort(diffs)[120] #Applying the Tricube Weight #...Tricube function tricube <- function(z) { ifelse (abs(z) < 1, (1 - (abs(z))^3)^3, 0) } #... a <- seq(0,1, by=.1) tricube(a) #Figure 2.5 plot(range(perotvote), c(0,1), xlab=\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, tricube((x.n - x0)/which.diff), cex=1) ########################################################################### #2、Kernel Smoothing ########################################################################### Figure 2.6 par(mfrow=c(3,1)) plot(perotvote, chal.vote, pch=\xlab=\main=\lines(ksmooth(perotvote, chal.vote, bandwidth=\plot(perotvote, chal.vote, pch=\xlab=\main=\lines(ksmooth(perotvote, chal.vote, kernel=\plot(perotvote, chal.vote, pch=\xlab=\main=\lines(ksmooth(perotvote, chal.vote, bandwidth=\#******* Kernel smoothing中选取box和normal核函数的比较，带宽相等 plot(perotvote, chal.vote, pch=\Share (%)\lines(ksmooth(perotvote, chal.vote, kernel=\lines(ksmooth(perotvote, chal.vote, kernel=\################################################################################## #第二部分，LPR模型 #Data Prep For Local Average Regression Step-by-Step cong <- as.data.frame(jacob[,2:3]) cong <- cong[order(cong$perotvote),1:2] y <- as.matrix(cong$chal.vote) x <- as.matrix(cong$perotvote) n <- length(y) #... tricube <- function(z) { ifelse (abs(z) < 1, (1 - (abs(z))^3)^3, 0) } #... x0 <- x[75] diffs <- abs(x - x0) which.diff <- sort(diffs)[120] x.n <- x[diffs<= which.diff] y.n <- y[diffs <= which.diff] weigh=tricube((x.n-x0)/which.diff) mod <- lm(y.n ~ x.n, weights=weigh) #Figure 2.7 plot(x, y, type=\bty=\abline(v=c(x0 - which.diff, x0 + which.diff), lty = 2) abline(v=x0) points(x[diffs > which.diff], y[diffs > which.diff], pch=16, cex=1, col=gray(.80)) points(x[diffs <= which.diff], y[diffs <= which.diff], cex=.85) abline(mod, lwd=2, col=1) text(27.5, 50, expression(paste(\x[0]))) #这里expression的用法比较有意思 arrows(25, 47, 15, 37, code =2, length = .10) ################################################################################# #2、Now Putting It Together For Local Regression Demonstration. #OLS Fit for Comparison ols <- lm(chal.vote ~ perotvote, data=jacob) #The loess fit model.loess <- loess(chal.vote ~ perotvote, data=jacob, span = 0.5) #*** 默认设置 degree=2，family=gauss, tricube加权 *** n <- length(chal.vote) x.loess <- seq(min(perotvote), max(perotvote), length=n) y.loess <- predict(model.loess, data.frame(perotvote=x.loess)) #得到预测值便于比较 #The lowess fit model.lowess <- lowess(chal.vote ~ perotvote, data=jacob, f = 0.5) #*** 默认设置 robust linear tricube加权 *** n <- length(chal.vote) x.lowess <- seq(min(perotvote), max(perotvote), length=n) y.lowess <- predict(model.lowess, data.frame(perotvote=x.lowess)) #得到预测值便于比较 #Figure 2.8 plot(perotvote, chal.vote, pch=\ylab=\lines(x.loess, y.loess) lines(x.lowess, y.lowess) abline(ols) legend(15,20, c(\################################################################################# #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 ~ perotvote, data=jacob, span = 0.5, degree=1, family=\iterations=1) #** no robust m1.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=\#*** 进行3次robust加权估计 *** plot(perotvote, chal.vote, pch=\(%)\lines(m1.lowess) lines(sort(perotvote), m1.loess$fit[order(perotvote)], lty=3, col=\lines(sort(perotvote), m0.loess$fit[order(perotvote)], lty=9,col=18) lines(m2.lowess, lty=2, col=\lines(sort(perotvote), m2.loess$fit[order(perotvote)], lty=4, col=\ ---------------------------------------------------------------------------- 第四章 样条估计spline 一、基本思想

按照x将样本分成多个区间，对每个区间分别进行估计。不同于核估计，这里不用移动计算，从而减小了计算量。 二、最简单的形式

Linear Spline with k knots:

其中，，

三、其他样条模型 1、p次样条估计

——二次样条Quadratic Spline (basis functions with k knots)

——三次样条Cubic Spline (with k knots, use quadratic basis functions)

——p-order spline (with k knots)

2、B-splines (with k knots cubic B-spline basis)

其中，

3、Natural Splines

以上估计方法对结点（knots）之间的估计比较准确，但对边界的拟合效果较差。自然样条的思想是，在自变量最小值和最大值处各增加一个结点，用线性模型拟合边界位置的样本点。

4、k的选择和模型比较 采用AIC准则

四、光滑样条smoothing spline 基于如果目标得到参数估计值

min

五、模型比较的F检验

六、R语言部分

library(foreign) jacob <- read.dta(\attach(jacob)

################################################# #第一部分，B样条和natural B 样条 library(splines) #*** P61 Perform Spline Regression **** m.bsp <- lm(chal.vote~bs(perotvote, df=5), data=jacob) #*** 3次B样条公式: df=k+3 (不含常数项) m.nsp <- lm(chal.vote~ns(perotvote, df=5), data=jacob) #*** df=5对应结点为4个；3次natural B样条公式：df=k+1 perot <- seq(min(perotvote), max(perotvote), length=312) bsfit <- predict(m.bsp, data.frame(perotvote=perot)) nsfit <- predict(m.nsp, data.frame(perotvote=perot)) AIC(m.bsp) #计算AIC值 ################################################################ #第二部分，光滑样条估计 plot(perotvote, chal.vote, pch=\ylab=\bty=\lines(smooth.spline(perotvote, chal.vote, df=2)) ################################################################## #第三部分，置信区间 library(splines) m.nsp <- lm(chal.vote~ns(perotvote, df=4), data=jacob) perot <- seq(min(perotvote), max(perotvote), length=312) nsfit <- predict(m.nsp, inteval=\ se.fit=TRUE, data.frame(perotvote=perot)) #Figure 3.8 plot(perotvote, chal.vote, pch=\(%)\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) #*** 偏移调整的置信区间 #*** Figure 3.9 #*** Overlay Natural cubic B-pline Fit and Confidence Intervals library(mgcv) m.smsp <- gam(chal.vote ~ s(perotvote, bs=\plot(m.smsp, rug=FALSE, se=TRUE, ylab=\residual=FALSE, shift=33.88) #绘制调整以后的上下界 lines(perot, nsfit$fit + 1.96*nsfit$se.fit, lty=3) lines(perot, nsfit$fit - 1.96*nsfit$se.fit, lty=3) legend(3,47, c(\lty=c(3, 2), bty=\################################################################## #第四部分，模型比较 ols <- lm(chal.vote ~ 1) #Use Automated R F-test Function Anova anova(ols, m.nsp) anova(ols, m.smsp) #** 三个不同的函数做smoothing spline，进行比较 plot(x,y, xlab=\#方法1 library(SemiPar) fit <- spm(y ~ f(x)) lines(fit, se=FALSE, lwd=1,lty=1,col=1) #方法2 library(pspline) lines(sm.spline(x,y, df=31),lty=5,col=5) #方法3 lines(smooth.spline(x,y,df=31),lty=6,col=6) legend(0,-0.8,c(\---------------------------------------------------------------------------- 第五章 Automated Smoothing Techniques自动光滑技术 一、Span by Cross-Validation

其中s指窗宽span。

【CV和GCV在LPR中表现不佳】 二、自动光滑技术

两种方法：1、采用MLE（极大似然估计）；2、采用CV选择。 1、MLE方法

得到。其中，是随机效应（the random effect）的方差估计，是随机误差

项（the error term）的方差估计 2、最小化CV或GCV

三、R语言部分

setwd(\课程/nonparameter regression/2010/ch3\jacob <- read.table(\attach(jacob) ########################################################################### #第一部分，计算GCV并绘图 library(locfit) alpha <- seq(0.2,0.8, by=0.05) plot(gcvplot(chal.vote~perotvote, data=jacob, alpha=alpha), type=\#第二部分，比较MLE方法和robust得到的估计结果 #***** P89 Figure 4.2 ***** plot(perotvote, chal.vote, pch=\ylab=\fit <- locfit(chal.vote~perotvote, data=jacob, alpha=0.2) lines(fit) lines(lowess(perotvote, chal.vote, f = 0.2), lty=2,col=2) # robust lines(lowess(perotvote, chal.vote, f = 0.2,iter=0), lty=3,col=3) # no robust legend(23,16, c(\lty=1:3,col=1:3,bty=\#第三部分，在光滑样条中是哟花姑娘GCV library(pspline) #使用其函数sm.spline() plot(perotvote, chal.vote, pch=\ylab=\main = \lines(sm.spline(perotvote, chal.vote, cv=F),lty=2,col=2) #GCV Smoothing Selection lines(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=\用LS方法估计，GCV smsp2 <- spm(chal.vote ~ f(perotvote)) # 光滑样条的混合模型表示，用极大似然方法估计！ smsp1 # or summary(smsp1) summary(smsp2) #***** P92 Figure 4.3 ***** par(mfrow=c(1,2)) plot(smsp1, rug=FALSE, se=FALSE, ylab=\residual=TRUE, shift=33.88, bty=\points(perotvote, chal.vote, pch=\plot(perotvote, chal.vote, pch=\ylab=\main=\lines(smsp2, se=FALSE, lwd=1, rug=FALSE)

---------------------------------------------------------------------------- 第六章 Additive and Semiparametric Regression Models 可加回归模型和半参数回归模型

R语言部分

jacob <- read.table(\attach(jacob) #第一部分，广义可加模型 library(mgcv) gam1 <- gam(chal.vote ~ s(perotvote, bs=\ + s(checks, bs=\, data=jacob) summary(gam1) #画图 par(mfrow = c(1,2)) plot(gam1, select=1, rug=F, se=TRUE, ylab=\bty=\points(perotvote, chal.vote, pch=\plot(gam1, select=2, rug=F, se=TRUE, ylab=\bty=\points(checks, chal.vote, pch=\#进行检查 gam.check(gam1) ##############--------------------------------------------------------------- #模型比较的卡方检验 #OLS Model ols1 <- gam(chal.vote ~ perotvote + checks, data=jacob) #或 ols <- lm(chal.vote ~ perotvote + checks, data=jacob) #速度更快 summary(gam1)$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) #第二种比较方法，采用anova anova(ols1, gam1, test=\########################################################## #第二部分，半参数模型 library(mgcv) #**** P123 Baseline Model ******** no transformations ****** ols <- gam(chal.vote ~ exp.chal + chal.spend + inc.spend + pres.vote + checks + marginal + partisan.redist + perotvote, data=jacob) #******* Test each continuous covariate gam1 <- gam(chal.vote ~ exp.chal + s(chal.spend, bs=\+ checks + marginal + partisan.redist + perotvote, data=jacob) gam2 <- gam(chal.vote ~ exp.chal + chal.spend + s(inc.spend, bs=\+ checks + marginal + partisan.redist + perotvote, data=jacob) gam3 <- gam(chal.vote ~ exp.chal + chal.spend + inc.spend + s(pres.vote, bs=\+ checks + marginal + partisan.redist + perotvote, data=jacob) gam4 <- gam(chal.vote ~ exp.chal + chal.spend + inc.spend + pres.vote + s(checks, bs=\+ 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=\anova(ols,gam1,test=\anova(ols,gam2,test=\anova(ols,gam3,test=\anova(ols,gam4,test=\anova(ols,gam5,test=\---------------------------------------------------------------------------- 第七章 Generalized Additive Models 广义可加模型

一、广义线性模型GLM

probit model:

;

logit model: ;

本文来源：https://www.bwwdw.com/article/3673.html