广义线性回归模型之泊松回归（logit/probit回归）—R语言实现

Poisson回归模型

Poisson回归也是广义线性回归模型中的一中，详细介绍可见之前的博客：
https://blog.csdn.net/qq_42871249/article/details/104339650

使用 Poisson回归模型时的两个问题

由于广义线性模型的假定很强, 所以当实际数据与假定的分布不符时会产生一些问题. Poisson回归模型也不例外, 人们目前主要关注的是 以下两个问题：

一、散布问题

在Poisson回归模型中, 假定方差和均值相等, 当方差大于或小于 均值时就会出现过散布(overdispersion)问题或欠散布 (underdispersion) 问题. 使用 Poisson回归模型时出现的散布问题的最简单解决办法是使用提到的准 Poisson回归模型, 而且还可以说明方差和均值的关系。准 Poisson 模型拟合代码示例：
glm(y~.,data,family=quasi(variance=“mu^2”,link=“log”))
这里的选项variance=“mu^2” 就把方差看成随着均值平方变化的函数, 这个选项可以输入"constant", “mu(1-mu)”, “mu”, “mu^2”, “mu^3”, 等 等.

二、零膨胀问题

代码介绍：

同样使用：glm(formula, family = gaussian, data, weights, subset, na.action, start = NULL, etastart, mustart, offset, control = glm.control(…), model = TRUE, method = “glm.fit”, x = FALSE, y = TRUE, contrasts = NULL, …)
不同之处在于参数family=possion()

仙客来案例分析

仙客来数据集包含三个自变量，因变量是Flowers开花数目。

直接上代码

w=read.csv("cyclamen.csv");n=nrow(w);m=ncol(w)
for(i in c(1,2,5))w[,i]=factor(w[,i])
table(w[,c(1,2,5)])

w=read.csv("cyclamen.csv");n=nrow(w)
for(i in c(1,2,5))w[,i]=factor(w[,i])
w=w[,-(3:4)]#不需要三行四行所以去掉
a=glm(Flowers~.,w,family=poisson)
summary(a)
anova(a,test="Chisq")#方差分析

b=step(a);summary(b)#用逐步回归法
anova(b,a,test="Chisq")

d=glm(b, w, family=quasipoisson(link = "log"))
anova(d,test="Chisq")

anova(d,test="Chisq")
anova(d,test="F")

glm(Flowers~.,data=w,family=quasi(variance="mu^2",link="log"))

下面做交叉验证：
首先定义了交叉验证的函数：Fold

Fold=function(Z=5,w,D,seed=7777){ 
  n=nrow(w);d=1:n;dd=list() 
  e=levels(w[,D]);T=length(e)#????��T??
  set.seed(seed)
  for(i in 1:T){
    d0=d[w[,D]==e[i]];j=length(d0)
    ZT=rep(1:Z,ceiling(j/Z))[1:j]
    id=cbind(sample(ZT,length(ZT)),d0);dd[[i]]=id}
 #上面每个dd[[i]]是随机1:Z及i类的下标集组成的矩阵
  mm=list()
  for(i in 1:Z){u=NULL;
  for(j in 1:T)u=c(u,dd[[j]][dd[[j]][,1]==i,2])
  mm[[i]]=u} #mm[[i]]Ϊ??i???±꼯i=1,...,Z
  return(mm)}#????Z???±꼯

w=read.csv("cyclamen.csv")
for(i in c(1,2,5))w[,i]=factor(w[,i])
w=w[,-(3:4)];D=4;Z=10
mm=Fold(Z,L,1,8888)#L[,1]??96??ˮƽ?ķ?????��

library(nnet)
MSE=matrix(99,Z,3);J=1
for(i in 1:Z)
{m=mm[[i]];M=mean((w[m,D]-mean(w[m,D]))^2)
a=glm(Flowers ~ ., family="poisson", data=w[-m,])
y1=predict(a,w[m,],type="response")
MSE[i,J]=mean((w[m,D]-y1)^2)/M}

J=J+1;for(i in 1:Z)
{m=mm[[i]];M=mean((w[m,D]-mean(w[m,D]))^2)
a=lm(Flowers ~ .,w[-m,])
MSE[i,J]=mean((w[m,D]-predict(a,w[m,]))^2)/M}

J=J+1;for(i in 1:Z)
  set.seed(1010);for(i in 1:Z){
    m=mm[[i]];M=mean((w[m,D]-mean(w[m,D]))^2)
    a=nnet(w[-m,D]/max(w[,D])~.,data=w[-m,],size=50,decay=0.1) 
    MSE[i,J]=mean((w[m,D]-predict(a,w[m,])*max(w[,D]))^2)/M}