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=\#*** 可以得到模型的每一项(terms)的预测值 sa2[,7] #*** f(nudem) 那一列 ##################################################################### #第二部分，用广义可加模型做多分类logit模型（顺序） library(foreign) couples <- read.dta(\课程/nonparameter regression/2010/Rprog/Chapter 6/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(parallel=T), data=couples) summary(gam1) par(mfrow=c(3,3)) plot(gam1,which.term=3, se=TRUE, rug=FALSE, bty=\plot(gam1,which.term=4, se=TRUE, rug=FALSE, bty=\plot(gam1,which.term=5, se=TRUE, rug=FALSE, bty=\plot(gam1,which.term=7, se=TRUE, rug=FALSE, bty=\plot(gam1,which.term=8, se=TRUE, rug=FALSE, bty=\plot(gam1,which.term=9, se=TRUE, rug=FALSE, bty=\plot(gam1,which.term=10, se=TRUE, rug=FALSE, bty=\【关于预测方法的说明没有看明白】 ##################################################################### #第三部分，用广义可加模型做possion模型 library(foreign) Scourt <- read.dta(\attach(Scourt) library(mgcv) m1 <- gam(nulls ~ tenure + congress + unified, data=Scourt,family=poisson) summary(m1) m2 <- gam(nulls ~ tenure + s(congress, bs=\family=poisson) summary(m2) anova(m1, m2, test=\#进行预测 predict.data <- data.frame(expand.grid(list(congress=seq(1, 101, by=1), tenure = 10.35, unified = 0))) #Created Fitted Values Using Fake Data 方法1 predict.fit <- predict(m2, newdata = predict.data, se.fit=TRUE) #Created Fitted Values Using Fake Data 方法2 predict.fit1 <- predict.gam(m1, newdata = predict.data, se.fit=TRUE) predict.fit2 <- predict.gam(m2, newdata = predict.data, se.fit=TRUE) #********************** Negative Binomial regression ********************** library(MASS) m4 <- gam(nulls~ tenure + s(congress, bs=\family=negative.binomial(1), control = gam.control(maxit = 150)) anova(m2, m4, test=\##################################################################### #第四部分，cox proportional hazards regression model【本部分的理论没有弄懂】 library(foreign) riot <- read.dta(\library(survival) m1 <- coxph(Surv(newend, cens) ~ lognwun + manuwage + unemrate + percfor + x15 + sumpi + lstwkall + lwka2 + d6768 + pasthist, data=riot) summary(m1) zph.m1 <- cox.zph(m1) m.base <- coxph(Surv(newend, cens) ~ pspline(nonwhtunemp, df=4) + pspline(manuwage, df=4) + pspline(unemrate, df=4) + pspline(percfor, df=4) + pspline(x15, df=4) + pspline(sumpi, df=4) + pspline(lstwkall, df=4) + d6768 + pasthist, data=riot) summary(m.base) #Test Nonproportional Hazards zph.m.base <- cox.zph(m.base) zph.m.base termplot(m.base, term=2, se=TRUE, rug=TRUE, ylab=\Wage\#Fit Parsimonious Model m.trim <- coxph(Surv(newend, cens) ~ pspline(nonwhtunemp, df=4) + pspline(manuwage, df=4) + unemrate + percfor + pspline(x15, df=4) + pspline(sumpi, df=4) + pspline(lstwkall, df=4) + d6768 + pasthist, data=riot) #Test Nonproportional Hazards zph.m.trim <- cox.zph(m.trim) zph.m.trim summary(m.trim) #Plot The Effects par(mfrow=c(2,3)) termplot(m.trim, term = 1, se = TRUE, ylab = \termplot(m.trim, term = 2, se = TRUE, ylab = \termplot(m.trim, term = 5, se = TRUE, ylab = \xlab = \termplot(m.trim, term = 6, se = TRUE, ylab = \termplot(m.trim, term = 7, se = TRUE, ylab = \#Figure 6.8 par(mfrow =c(1,2)) termplot(m.trim, term = 2, se = TRUE, ylab = \termplot(m.trim, term = 1, se = TRUE, ylab = \anova(m.trim, m.base, test=\