################################################# #第一部分，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: ;