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

三、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: ;

possion regression:

二、广义可加模型

或

三、估计方法

MLE: use Newton-Raphson algorithm IRLS: backfitting algorithm (in ch5) 四、假设检验

LR=-2（LogLikelihood0 – LogLikelihood1） 【这是两个模型进行比较】 五、R语言部分

setwd() library(foreign) war <- read.dta(\attach(war) #第一部分，用广义可加模型做Logit回归【与glm的用法区别】 library(mgcv) war.glm <- gam(dispute ~ nudem + nugrow + allies + contig + nucapab + trade + sumdisp + s(year, bs=\family=binomial) summary(war.glm) #进行预测 war.smooth7 <- gam(dispute ~ s(nudem, bs=\+ sumdisp + s(year, bs=\summary(war.smooth7) new <- data.frame(nudem=nudem, nugrow=nugrow, allies=allies, contig=contig, 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=\