######################  导入所需程序包
library(data.table)
library(ggplot2)
library(hglm)
library(car)
library(lmtest)


######################  导入数据
data=read.csv("D:\\Other\\Personal\\研究生案例大赛\\研究生案例大赛\\collision_data.csv",header=T)
data=data.table(data)
names(data)


#####################  单因素分析
data[,sum(车损赔款),by='省份']
data[,sum(车损赔款),by='车龄']
data[,sum(车损赔款),by='载客数']
min(data$车损赔款);max(data$车损赔款);mean(data$车损赔款)
qplot(data$车损赔款,geom="histogram",xlim=c(0,10000),xlab="车损赔款")
ss=sort(data$车损赔款)
n=length(ss)
length(ss[ss<=3000])/n


#####################  广义线性模型GLM
##  导入数据
data=read.csv("e:\\数据\\collision_data.csv",header=T)
data=data.table(data)
names(data)

##  选择基准水平
data$省份=relevel(factor(data$省份),ref='山东')#
data$车龄=relevel(factor(data$车龄),ref="0.5")
data$载客数=relevel(factor(data$载客数),ref="5")

##  广义线性模型
glm1=glm(车损赔款~省份+车龄+载客数,data=data,family=Gamma(link="log"))
Anova(glm1,type="III")
summary(glm1)

glm2=glm(车损赔款~省份,data=data,family=Gamma(link="log"))
Anova(glm2,type="III")
summary(glm2)

##  计算拟合值
XX=model.matrix(~省份,data=data)
yy=data$车损赔款
fitted_glm=exp(XX%*%(glm2$coef))
省份=levels(data$省份)
车龄=levels(data$车龄)
载客数=levels(data$载客数)
fit1=rep(0,length(省份))
fit2=rep(0,length(车龄))
fit3=rep(0,length(载客数))
names(fit1)=省份
names(fit2)=车龄
names(fit3)=载客数
for (i in 1:length(省份)){
  fit1[i]=mean(fitted_glm[data$省份==省份[i]])
}
for(i in 1:length(车龄)){
  fit2[i]=mean(fitted_glm[data$车龄==车龄[i]])
}
for(i in 1:length(载客数)){
  fit3[i]=mean(fitted_glm[data$载客数==载客数[i]])
}
glm_fit=c(fit1,fit2,fit3)


##  计算实际值
省份=levels(data$省份)
车龄=levels(data$车龄)
载客数=levels(data$载客数)
act1=rep(0,length(省份))
act2=rep(0,length(车龄))
act3=rep(0,length(载客数))
names(act1)=省份
names(act2)=车龄
names(act3)=载客数

for (i in 1:length(省份)){
  act1[i]=mean(data$车损赔款[data$省份==省份[i]])
}
for(i in 1:length(车龄)){
  act2[i]=mean(data$车损赔款[data$车龄==车龄[i]])
}
for(i in 1:length(载客数)){
  act3[i]=mean(data$车损赔款[data$载客数==载客数[i]])
}
act=c(act1,act2,act3)


#拟合值与实际值的对比
cbind(glm_fit,act)



######################  广义线性混合模型GLMM
##  导入数据
data=read.csv("e:\\数据\\collision_data.csv",header=T)
data=data.table(data)
names(data)

##  选择基准水平
data$省份=relevel(factor(data$省份),ref='山东')#
data$车龄=relevel(factor(data$车龄),ref="0.5")
data$载客数=relevel(factor(data$载客数),ref="5")

##  随机效应独立，且服从正态分布
XX=model.matrix(~车龄+载客数,data=data)
yy=data$车损赔款
ZZ=model.matrix(~factor(省份)-1,data=data)
glmm=hglm(X=XX,y=yy,Z=ZZ,family=Gamma(link=log),
           rand.family=gaussian(link="identity"),data=data,calc.like=T)
summary(glmm)
lrt(glmm)


##  计算拟合值
fitted_glmm=exp(XX%*%(glmm$fixef)+ZZ%*%(glmm$ranef))
省份=levels(data$省份)
车龄=levels(data$车龄)
载客数=levels(data$载客数)
fit1=rep(0,length(省份))
fit2=rep(0,length(车龄))
fit3=rep(0,length(载客数))
names(fit1)=省份
names(fit2)=车龄
names(fit3)=载客数
for (i in 1:length(省份)){
  fit1[i]=mean(fitted_glmm[data$省份==省份[i]])
}
for(i in 1:length(车龄)){
  fit2[i]=mean(fitted_glmm[data$车龄==车龄[i]])
}
for(i in 1:length(载客数)){
  fit3[i]=mean(fitted_glmm[data$载客数==载客数[i]])
}
glmm_fit=c(fit1,fit2,fit3)


##  计算实际值
省份=levels(data$省份)
车龄=levels(data$车龄)
载客数=levels(data$载客数)
act1=rep(0,length(省份))
act2=rep(0,length(车龄))
act3=rep(0,length(载客数))
names(act1)=省份
names(act2)=车龄
names(act3)=载客数

for (i in 1:length(省份)){
  act1[i]=mean(data$车损赔款[data$省份==省份[i]])
}
for(i in 1:length(车龄)){
  act2[i]=mean(data$车损赔款[data$车龄==车龄[i]])
}
for(i in 1:length(载客数)){
  act3[i]=mean(data$车损赔款[data$载客数==载客数[i]])
}
act=c(act1,act2,act3)


#拟合值与实际值的对比
cbind(glmm_fit,act)




###############  分层广义线性模型hglm
###  导入数据
data=read.csv("e:\\数据\\collision_data.csv",header=T)
data=data.table(data)
names(data)


###  选择基准水平
data$省份=relevel(factor(data$省份),ref='山东')#
data$车龄=relevel(factor(data$车龄),ref="0.5")
data$载客数=relevel(factor(data$载客数),ref="5")


#####分层广义线性模型：随机效应独立，服从伽马分布
XX=model.matrix(~车龄+载客数,data=data)
yy=data$车损赔款
ZZ=model.matrix(~factor(省份)-1,data=data)
hglm=hglm(X=XX,y=yy,Z=ZZ,family=Gamma(link="log"),
           rand.family=Gamma(link="identity"),data=data,calc.like=T)
summary(hglm)
lrt(hglm)


##  计算拟合值
fitted_hglm=exp(XX%*%(hglm$fixef)+ZZ%*%(hglm$ranef))
省份=levels(data$省份)
车龄=levels(data$车龄)
载客数=levels(data$载客数)
fit1=rep(0,length(省份))
fit2=rep(0,length(车龄))
fit3=rep(0,length(载客数))
names(fit1)=省份
names(fit2)=车龄
names(fit3)=载客数

for (i in 1:length(省份)){
  fit1[i]=mean(fitted_hglm[data$省份==省份[i]])
}
for(i in 1:length(车龄)){
  fit2[i]=mean(fitted_hglm[data$车龄==车龄[i]])
}
for(i in 1:length(载客数)){
  fit3[i]=mean(fitted_hglm[data$载客数==载客数[i]])
}
hglm_fit=c(fit1,fit2,fit3)


#拟合值与实际值的对比
cbind(hglm_fit,act)



#################  空间分层广义线性模型shglm
##  导入数据
data=read.csv("e:\\数据\\collision_data.csv",header=T)
data=data.table(data)
names(data)


##  选择基准水平
data$省份=relevel(factor(data$省份),ref='山东')#
data$车龄=relevel(factor(data$车龄),ref="0.5")
data$载客数=relevel(factor(data$载客数),ref="5")


##  随机效应服从联合正态分布
XX=model.matrix(~车龄+载客数,data=data)
yy=data$车损赔款
ZZ=model.matrix(~factor(省份)-1,data=data)
DD=matrix(c(rep(0,13),rep(1,2),rep(0,2),1,0,1,rep(0,2),rep(1,2),0),
          nrow=5,byrow=T)
DD
Rnd=ncol(ZZ)
shglm=hglm(X=XX,y=yy,Z=ZZ,family=Gamma(link="log"),
                     rand.family=CAR(D=DD),data=data,calc.like=T)
summary(shglm)
lrt(shglm)


##  计算拟合值
fitted_shglm=exp(XX%*%(shglm$fixef)+ZZ%*%(shglm$ranef))
省份=levels(data$省份)
车龄=levels(data$车龄)
载客数=levels(data$载客数)
fit1=rep(0,length(省份))
fit2=rep(0,length(车龄))
fit3=rep(0,length(载客数))
names(fit1)=省份
names(fit2)=车龄
names(fit3)=载客数

for (i in 1:length(省份)){
  fit1[i]=mean(fitted_shglm[data$省份==省份[i]])
}
for(i in 1:length(车龄)){
  fit2[i]=mean(fitted_shglm[data$车龄==车龄[i]])
}
for(i in 1:length(载客数)){
  fit3[i]=mean(fitted_shglm[data$载客数==载客数[i]])
}
shglm_fit=c(fit1,fit2,fit3)


##  拟合值与实际值对比
cbind(shglm_fit,act)




######## 模型对比
exposure=c(table(data$省份),table(data$车龄),table(data$载客数))
cbind(exposure,glmm_fit,hglm_fit,shglm_fit,act)
cbind(glmm$likelihood[1],hglm$likelihood[1],shglm$likelihood[1])