R语言学习笔记

文章目录

R语言学习笔记
probit回归
- factor()和as.factor()
- relevel()
案例11.4复刻 glm函数
- 整理变量
- 回归：Logistic和Probit——glm()
- - family=binomial(link=logit)
  - family=binomial(link=probit)
- 计算log likelihood
- 结果输出
比较probit和logistic

R语言中的数据类型之factor因子

probit回归

factor()和as.factor()

这两个没区别。

#MEPS DATA
Hexpend<-read.csv("HealthExpend.csv")  #导入数据# CHECK THE NAMES,DIMENSION IN THE FILE AND LIST THE FRIST
names(Hexpend)
dim(Hexpend)
Hexpend[1:8,]
attach(Hexpend)
n<-dim(Hexpend)[1]
POSEXP<-seq(0,0,length=n)
for(i in 1:n){if(EXPENDIP[i]!=0)POSEXP[i]=1}# ALTERNATIVE - FIT A GENERALIZED LINEAR MODEL;
PosExpglm = glm(POSEXP~GENDER,family=binomial(link=logit))
summary(PosExpglm)
logLik(PosExpglm)
summary(POSEXP)# FULL LOGIT MODEL
PosExpglmFull=glm(POSEXP~AGE+GENDER+factor(RACE)+factor(REGION)+factor(EDUC)+factor(PHSTAT)+factor(ANYLIMIT)+factor(INCOME)+factor(insure),family=binomial(link=logit))
summary(PosExpglmFull)
logLik(PosExpglmFull)Gender<-as.factor(GENDER)
PosExpglmFull=glm(POSEXP~AGE+C(Gender,base=1)+as.factor(RACE)+as.factor(REGION)+as.factor(EDUC)+as.factor(PHSTAT)+as.factor(ANYLIMIT)+as.factor(INCOME)+as.factor(insure),family=binomial(link=logit))
summary(PosExpglmFull)

relevel()

# CHANGE REFERANCE LEVELS TO AGREE WITH BOOK(DONE IN SAS)
RACE=relevel(factor(RACE),ref="WHITE")
REGION=relevel(factor(REGION),ref="WEST")
EDUC=relevel(factor(EDUC),ref="LHIGHSC")
PHSTAT=relevel(factor(PHSTAT),ref="EXCE")
INCOME=relevel(factor(INCOME),ref="POOR")

这里给的例子里没有factor，但不加会报错。
可以对比relevel前后的结果

案例11.4复刻 glm函数

Regression Modeling with Actuarial and Financial P315 案例11.4 Application: Medical Expenditures

下载数据HealthExpend.csv

HealthExpend<-read.csv("HealthExpend.csv")  #导入数据
str(HealthExpend)attach(HealthExpend)
n=2000 #共有2000个数据sink("result.txt")  #回归结果等储存在这里

整理变量

## 整理变量
#Ethnicity即RACE或RACE1
ASIAN<-seq(0,0,length=n)    # 1 if Asian 4.3 4.7
BLACK<-seq(0,0,length=n)    # 1 if Black 14.8 10.5
NATIVE<-seq(0,0,length=n)   #1 if Native 1.1 13.6???这个的描述统计结果不符#WHITE Reference level 79.9 7.5???other怎么处理#Region即REGION或REGION1
NORTHEAST<-seq(0,0,length=n)    # 1 if Northeast 14.3 10.1
MIDWEST<-seq(0,0,length=n)  # 1 if Midwest 19.7 8.7
SOUTH<-seq(0,0,length=n)    #1 if South 38.2 8.4#WEST Reference level 27.9 5.4#Education 即EDUC
COLLEGE<-seq(0,0,length=n)  # 1 if college or higher degree 27.2 6.8
HIGHSCHOOL<-seq(0,0,length=n)   # 1 if high school degree 43.3 7.9#Reference level is lower than high school degree 29.5 8.8#Self-rated physical health 即PHSTAT和PHSTAT1
POOR<-seq(0,0,length=n) # 1 if poor 3.8 36.0
FAIR<-seq(0,0,length=n) # 1 if fair 9.9 8.1
GOOD<-seq(0,0,length=n) # 1 if good 29.9 8.2
VGOOD<-seq(0,0,length=n)    # 1 if very good 31.1 6.3#Reference level is excellent health 25.4 5.1#Income compared to poverty line 即INCOME和INCOME1
HINCOME<-seq(0,0,length=n)  # 1 if high income 31.6 5.4
MINCOME<-seq(0,0,length=n)  # 1 if middle income 29.9 7.0
LINCOME<-seq(0,0,length=n)  # 1 if low income 15.8 8.3
NPOOR<-seq(0,0,length=n)    # 1 if near poor 5.8 9.5#Reference level is poor/negative 17.0 13.0#Insurance coverage insurance in any month of 2003
INSURE=insure# 1 if covered by public/private health 77.8 9.2 0 if have no health insurance in 2003 22.3 3.1#因变量
Y<-seq(0,0,length=n)for(i in 1:n)
{if(RACE[i]=="ASIAN")   ASIAN[i]=1if(RACE[i]=="BLACK") BLACK[i]=1if(RACE[i]=="NATIV") NATIVE[i]=1if(REGION[i]=="NORTHEAST")  NORTHEAST[i]=1if(REGION[i]=="MIDWEST") MIDWEST[i]=1   if(REGION[i]=="SOUTH")  SOUTH[i]=1 if(EDUC[i]=="COLLEGE")  COLLEGE[i]=1   if(EDUC[i]=="HIGHSCH")  HIGHSCHOOL[i]=1if(PHSTAT[i]=="POOR")   POOR[i]=1  if(PHSTAT[i]=="FAIR")   FAIR[i]=1if(PHSTAT[i]=="GOOD") GOOD[i]=1if(PHSTAT[i]=="VGOO") VGOOD[i]=1if(INCOME[i]=="HINCOME") HINCOME[i]=1if(INCOME[i]=="MINCOME")   MINCOME[i]=1if(INCOME[i]=="LINCOME")   LINCOME[i]=1if(INCOME[i]=="NPOOR") NPOOR[i]=1if(EXPENDIP[i]!=0)  Y[i]=1
}data<-data.frame(Y,AGE,GENDER,ASIAN,BLACK,NATIVE,NORTHEAST,MIDWEST,SOUTH,COLLEGE,HIGHSCHOOL,POOR,FAIR,GOOD,VGOOD,MNHPOOR,ANYLIMIT,HINCOME,MINCOME,LINCOME,NPOOR,INSURE) #用于回归的数据集
detach(HealthExpend)
attach(data)
m<-length(data)  #变量数
summary(data)   #描述统计

回归：Logistic和Probit——glm()

family=binomial(link=logit)

Logistic.full<-glm(Y~.,data=data,family=binomial(link=logit))
print(summary(Logistic.full))Logistic.reduced<-glm(Y~.-AGE-MNHPOOR,data=data,family=binomial(link=logit))
print(summary(Logistic.reduced))

family=binomial(link=probit)

Probit.reduced<-glm(Y~.-AGE-MNHPOOR,data=data,family=binomial(link=probit))
print(summary(Probit.reduced))

计算log likelihood

这个可以直接用logLik函数计算，这里是自己手动写了个程序……
我也是对自己无语了哈哈哈

#计算log likelihood===============================
LogLikelihood<-c(0,0,0)
logit<-function(z){ 1/(1+exp(-z))
}
#Logistic.full-----------------------------------
result<-seq(0,0,length=2000)
z<-seq(0,0,length=2000)
pai<-seq(0,0,length=2000)for(i in 1:n){for(j in 1:(m-1)){z[i]<-coef(Logistic.full)[j+1]*data[i,j+1]+z[i]}z[i]=z[i]+coef(Logistic.full)[1]pai[i]=logit(z[i])result[i]=Y[i]*log(pai[i])+(1-Y[i])*log(1-pai[i])
}LogLikelihood[1]<-sum(result)#Logistic.reduced-----------------------------------
result<-seq(0,0,length=2000)
z<-seq(0,0,length=2000)
pai<-seq(0,0,length=2000)for(i in 1:n){for(j in 1:(m-3)){z[i]<-coef(Logistic.reduced)[j+1]*data[i,names(coef(Logistic.reduced)[j+1])]+z[i]}z[i]=z[i]+coef(Logistic.reduced)[1]pai[i]=logit(z[i])result[i]=Y[i]*log(pai[i])+(1-Y[i])*log(1-pai[i])
}LogLikelihood[2]<-sum(result)#Probit.reduced-----------------------------------
result<-seq(0,0,length=2000)
z<-seq(0,0,length=2000)
pai<-seq(0,0,length=2000)for(i in 1:n){for(j in 1:(m-3)){z[i]<-coef(Probit.reduced)[j+1]*data[i,names(coef(Probit.reduced)[j+1])]+z[i]}z[i]=z[i]+coef(Probit.reduced)[1]pai[i]=pnorm(z[i])result[i]=Y[i]*log(pai[i])+(1-Y[i])*log(1-pai[i])
}
LogLikelihood[3]<-sum(result)cat("LogLikelihood: ",LogLikelihood)sink()

结果输出

#系数
write.csv(summary(Logistic.full)$coef,file="Logistic_full.csv")
write.csv(summary(Logistic.reduced)$coef,file="Logistic_reduced.csv")
write.csv(summary(Probit.reduced)$coef,file="Probit_reduced.csv")#描述统计
col.name<-c("min","max","median","mean","sd")
describe<-matrix(data = NA, nrow = m, ncol = 5, byrow = FALSE,dimnames =list(names(data),col.name))for(i in 1:m){describe[i,1]=min(data[,i])describe[i,2]=max(data[,i])describe[i,3]=median(data[,i])describe[i,4]=mean(data[,i])describe[i,5]=sd(data[,i])
}write.csv(data.frame(describe),file="summary.csv")

比较probit和logistic

#Probit模型
curve(exp(x)/(1+exp(x)),xlim=c(-4,4),col="red",ylab="Distribution Function")
text(2.3,.75,"Logit Case",cex=1.2,col="red")
arrows(2.3,.79,1.88,.84,length=0.1,angle=20,col="red")#Logit模型
curve(pnorm(x),xlim=c(-4,4),col="blue",add=T,ylab="Distribution Function")
text(0,.95,"Probit Case",cex=1.2,col="blue")
arrows(-.1,0.9,.6,.8,length=0.1,angle=20,col="blue")

R语言学习笔记 07 Probit、Logistic回归相关推荐

R语言学习笔记 06 岭回归、lasso回归
R语言学习笔记文章目录 R语言学习笔记比较lm.ridge和glmnet函数画岭迹图图6-4 <统计学习导论基于R语言的应用>P182 图6-6<统计学习导论基于R语言的 ...
R语言学习笔记——高级篇：第十四章-主成分分析和因子分析
R语言 R语言学习笔记--高级篇:第十四章-主成分分析和因子分析文章目录 R语言前言一.R中的主成分和因子分析二.主成分分析 2.1.判断主成分的个数 2.2.提取主成分 2.3.主成分旋转 ...
R语言splines包构建基于logistic回归的自然样条分析：南非心脏病数据集、非线性：基函数展开和样条分析、你简单分析的不重要特征，可能只是线性不显著、而非线性是显著的
R语言splines包构建基于logistic回归的自然样条分析:南非心脏病数据集.非线性:基函数展开和样条分析.你简单分析的不重要特征,可能只是线性不显著.而非线性是显著的目录
R语言学习笔记（1~3）
R语言学习笔记(1~3) 一.R语言介绍 x <- rnorm(5) 创建了一个名为x的向量对象,它包含5个来自标准正态分布的随机偏差. 1.1 注释由符号#开头. #函数c()以向量的形式输 ...
r语言c函数怎么用,R语言学习笔记——C#中如何使用R语言setwd()函数
在R语言编译器中,设置当前工作文件夹可以用setwd()函数. > setwd("e://桌面//") > setwd("e:\桌面\") > ...
R语言学习笔记——入门篇：第一章-R语言介绍
R语言 R语言学习笔记--入门篇:第一章-R语言介绍文章目录 R语言一.R语言简介 1.1.R语言的应用方向 1.2.R语言的特点二.R软件的安装 2.1.Windows/Mac 2.2.Lin ...
R语言学习笔记——入门篇：第三章-图形初阶
R语言 R语言学习笔记--入门篇:第三章-图形初阶文章目录 R语言一.使用图形 1.1.基础绘图函数:plot( ) 1.2.图形控制函数:dev( ) 补充--直方图函数:hist( ) 补充- ...
R语言学习笔记（八）--读写文件与网络爬虫
R语言学习笔记(八) 1 工作路径 2 保存R对象 3 Scan函数 3-1 从控制台读取数据 3-2 从txt文件读取数据 3-3 从url读取数据 4 按行读写文本文件 5 读取文本文件(txt. ...
R语言学习笔记（三）多元数据的数据特征、相关分析与图形表示
文章目录写在前面独立性检验 χ2\chi^2χ2独立性检验 Fisher独立性检验 Cochran-Mantel-Haenszel χ2\chi^2χ2独立性检验相关性分析相关性检验相关性检 ...

R语言学习笔记 07 Probit、Logistic回归