新人求助，如何使用R软件模拟probit模型数据？悬赏1500币，满意可以追加。

本人想构建一个双变量Probit模型来估计二分类数据的回归系数。看着很长，问题很简单。因为数据具有内生性，所以会用到工具变量，构建两个Probit模型的回归方程。现在在数据模拟时遇到问题。
   思路比较简单，不知道对不对。
   程序中g1为工具变量，c1、c2为两个外生变量，u为模拟的遗漏的混杂因素，在数据模拟时会用到，后续参数估计时会剔除来
   估计模型的效果，exi为残差项。
    g1=rbinom(10000, 1, 0.9)                    #模拟数据
    c1=rnorm(10000, mean=10, sd=2)
    c2=rbinom(10000, 2, 0.4)
    u=rnorm(10000, mean = 10, sd = 2)
    exi=rnorm(10000, mean = 0, sd = 1)
    eyi=rnorm(10000, mean = 0, sd = 1)
    首先、根据自变量和设定的参数计算出内生变量X,不知道这样确定的X是否正确。用R的glm模型估计出来的参数与模拟参数很接近。
    probit_x=0.3*g1+0.2*c1-0.3*c2+u+exi-13
    probit_data=as.data.frame(cbind(g1,c1,c2,u,probit_x))
    probit_data$x[probit_data$probit_x>0]=1
    probit_data$x[probit_data$probit_x<=0]=0
    anorex.1 <- glm(x~g1+c1 +c2+u,
                    family=binomial(link = "probit"),data =probit_data)
    summary(anorex.1)
    第二步、根据第一步模拟的内生变量X和外生变量，计算出因变量Y.但是这一次估计的系数与模拟系数差异太大。
    probit_y=0.3*x+2*c1-4*c2+*u-13
    probit_data=as.data.frame(cbind(probit_data,probit_y))
    probit_data$y[probit_data$probit_y>0]=1
    probit_data$y[probit_data$probit_y<=0]=0
    anorex.2 <- glm(y~x+c1 +c2+u,
                    family=binomial(link = "probit"),data =probit_data)
    summary(anorex.2)

    研究的最终目的是估计内生变量的参数与模拟的是否一样。求各位大神指出数据模拟时的问题。


    此外，后续将在R软件中用贝叶斯方法      估计在去掉遗漏因素U（构建方程中不会带入U）的情况下，内生变量X的回归参数
与模拟的差别有多大。模型为probit模型，以及双变量Probit模型。如果有高手，请发邮件至lnxiangchun@163.com，或
站内信，有酬劳，希望数学系或计量经济学大牛救急，毕业课题做不下去了啊。

1. #第二部分代码修改如下
   
2. x<-probit_x
   
3. probit_y=0.3*x+2*c1-4*c2+u+eyi-13
   
4.     probit_data=as.data.frame(cbind(probit_data,probit_y))
   
5.     probit_data$y[probit_data$probit_y>0]=1
   
6.     probit_data$y[probit_data$probit_y<=0]=0
   
7.     anorex.2 <- glm(y~x+c1 +c2+u,
   
8.                     family=binomial(link = "probit"),data =probit_data)
   
9.     summary(anorex.2)

复制代码

统计R浪人 发表于 2013-10-12 21:17

谢谢你的回答，不过在我的模拟过程中，X和Y均为二分类变量，是通过设定的标准参数以及自变量通过probit模型计算出来的，后续分析估计将剔除未知混杂因素U来计算回归参数，看回归参数与原来模拟的标准参数的差异。

没有人么，自己顶。求助 如何在R进行probit模型数据的模拟？

如何根据设定的回归参数及自变量模拟金标准因变量。

楼主就7个论坛币，您是怎么悬赏1500论坛币的？想骗答案啊？

data <- read.csv("I:/R/R2/kechengsheji.csv")    #读取数据
a1 <- data[1:50,]          #取前50行的数据建立模型
a2 <- data[51:60,]         #取后10行的数据进行预测
par(mfrow = c(2,2))        #设置画图格式为2*2
boxplot(sex~yonxhsx,data = a1,main ="sex")
boxplot(zy~yonxhsx,data = a1,main ="zy")
boxplot(ah~yonxhsx,data = a1,main ="ah")
boxplot(swsj~yonxhsx,data = a1,main ="swsj")
par(mfrow = c(2,2))
boxplot(yxsj~yonxhsx,data = a1,main ="yxsj")
boxplot(zyqk~yonxhsx,data = a1,main ="zyqk")
boxplot(sxcj~yonxhsx,data = a1,main ="sxcj")
glm0.a <- glm(yonxhsx~1,family = binomial(link = logit),data = a1)       #拟合logit回归，不使用任何变量的空模型
glm1.a <- glm(yonxhsx~sex+zy+ah+swsj+yxsj+zyqk+sxcj, family = binomial(link = logit),data = a1)  #logit回归全模型
anova(glm0.a,glm1.a)     #计算glm0.a和glm1.a的dexiance
1-pchisq(33.978,6)       #计算模型的显著性检验的P值
glm0.b <- glm(yonxhsx~1,family = binomial(link = probit),data = a1)      #拟合probit回归，不使用任何变量的空模型
glm1.b <- glm(yonxhsx~sex+zy+ah+swsj+yxsj+zyqk+sxcj, family = binomial(link = probit),data = a1)    #probit回归全模型
anova(glm0.b,glm1.b)    #计算glm0.b和glm1.b的dexiance
1-pchisq(33.834,6)      #计算模型的显著性检验的P值
library(car)            #载入程序包car
anova(glm1.a,type="lll")   #对模型glm1.a作三型方差分析
summary(glm1.a)   #显示模型glm1.a的各方面细节，包括参数估计值、P值等
anova(glm1.b,type="lll")    #对模型glm1.b作三型方差分析
summary(glm1.b)    #显示模型glm1.b的各方面细节，包括参数估计值、P值等
glm2.b <- glm(yonxhsx~zy+ah+zyqk, family = binomial(link = probit),data = a1)     #拟合简化后的probit 回归模型
summary(glm2.b)      #显示模型glm2.b的各方面细节，包括参数估计值、P值等
glm2.a <- glm(yonxhsx~zy+ah+zyqk, family = binomial(link = logit),data = a1)    #拟合简化后的logit 回归模型
summary(glm2.a)    #显示模型glm2.a的各方面细节，包括参数估计值、P值等
deviance(glm0.a)       #计算模型glm0.a的deviance
deviance(glm1.a)       #计算模型glm1.a的deviance
AIC(glm0.a)            #计算模型glm0.a的AIC取值
AIC(glm1.a)            #计算模型glm1.a的AIC取值
AIC(glm0.a,k = log(length(a1[,1])))       #计算模型glm0.a的BIC取值
AIC(glm1.a,k = log(length(a1[,1])))       #计算模型glm1.a的BIC取值
logit.aic <- step(glm1.a,trace=0)    #根据AIC准则选择最优模型，并赋值给logit.aic
summary(logit.aic)    #显示模型logit.aic的各方面细节，包括参数估计值、P值等
n <- length(a1[,1])    #样本大小
logit.bic <- step(glm1.a,k = log(n),trace=0)    #根据BIC准则选择最优模型，并赋值给logit.bic
summary(logit.bic)       #显示模型logit.bic的各方面细节，包括参数估计值、P值等
probit.aic = step(glm1.b,trace = 0)    #根据AIC准则选择最优模型，并赋值给probit.aic
summary(probit.aic)     #显示模型probit.aic的各方面细节，包括参数估计值、P值等
probit.bic = step(glm1.b,k = log(n),trace = 0)   #根据BIC准则选择最优模型，并赋值给probit.bic
summary(probit.bic)     #显示模型probit.bic的各方面细节，包括参数估计值、P值等
summary(glm1.a)     #显示模型glm1.a的各方面细节，包括参数估计值、P值等
p <- predict(glm2.a,a2)    #利用模型glm1.a对数据a2进行预测
p <- exp(p)/(1+exp(p))    #计算预测得到的概率
a2$yonxhsx.pred <- 1*(p>0.5)   #以0.5为阈值生成预测值
table(a2[,c(7,8)])           #计算预测值与真实值的2维数表
a2$yonxhsx.pred <- 1*(p>0)    #以0为阈值生成预测值
table(a2[,c(7,8)])           #计算预测值与真实值的2维数表
a2$yonxhsx.pred <- 1*(p>0.05)   #以0.05为阈值生成预测值
table(a2[,c(7,8)])         #计算预测值与真实值的2维数表
ngrids <- 100           #设置格点数为100
TPR <- rep(0,ngrids)    #为TPR（true positive rate）赋初值
FPR <- rep(0,ngrids)    #为FPR（false positive rate）赋初值
for(i in 1:ngrids){
        p0 = i/ngrids                #选取初值p0
        yonxhsx.true = a2$yonxhsx    #从a2中选出真实值并赋给yonxhsx.true
        yonxhsx.pred = 1*(p > p0)      #以0.05为阈值生成预测值
        FPR[i] = sum(yonxhsx.pred*(1-yonxhsx.true))/sum(1-yonxhsx.true)   #计算FPR
        TPR[i] = sum(yonxhsx.pred*yonxhsx.true)/sum(yonxhsx.true)  #计算TPR
}
plot(FPR,TPR,type = "l",col = 2)    画出FPR与TPR的散点图，即ROC曲线
points(c(0,1),c(0,1),type = "l",lty = 2)    #添加对角线
p <- matrix(0,length(a2[,1]),6)    #生成矩阵，用于存储各模型的预测值
p[,1] <- predict(glm1.a,a2)      #利用模型glm1.a对数据a2进行预测
p[,2] <- predict(logit.aic,a2)    #利用模型logit.aic对数据a2进行预测
p[,3] <- predict(logit.bic,a2)    #利用模型logit.bic对数据a2进行预测
p[,c(1:3)] <- exp(p[,c(1:3)])/(1+exp(p[,c(1:3)]))  #计算所得到的概率
p[,4] <- predict(glm1.b,a2)   #利用模型glm1.b对数据a2进行预测
p[,5] <- predict(logit.aic,a2)   #利用模型logit.aic对数据a2进行预测
p[,6] <- predict(logit.bic,a2)   #利用模型logit.bic对数据a2进行预测
p[,c(4:6)] <- exp(p[,c(4:6)])/(1+exp(p[,c(4:6)]))    #计算所得到的概率
plot(c(0,1),c(0,1),type = "l",main = "FPR vs TPR",xlab = "FPR",ylab = "TPR")      #画图，生成基本框架
TPR <- rep(0,ngrids)     #为TPR赋初值
FPR <- rep(0,ngrids)     #为FPR赋初值
for(k in 1:6){
        prob = p[,k]  #取出p中的k列的值，即第k个模型的预测概率
        for(i in 1:ngrids){
                p0 = i/ngrids         #选取阈值
                yonxhsx.hat = 1*(prob > p0)    #根据阈值生成预测值
                TPR[i] = sum(yonxhsx.true*yonxhsx.hat)/sum(yonxhsx.true)   #计算TPR
                FPR[i] = sum(yonxhsx.hat*(1-yonxhsx.true))/sum(1-yonxhsx.true)   #计算FPR
        }
        points(FPR,TPR,type = "b",col = k,lty = k,pch = k)   #向图上添加第k个模型的TPR与FPR的散点图
}
legend(0.6,0.5,c("LOGIT FULL MODEL","LOGIT AIC MODEL",
"LOGIT BIC MODEL","PROBIT FULL MODEL","PROBIT AIC MODEL",
"PROBIT BIC MODEL"),lty = c(1:6),col = c(1:6),pch = c(1:6)) #为6个模型添加标识，区分6个模型

好久了。这个问题好像还是没有解决。。。谢谢各位啦。

烟草味99 发表于 2020-3-27 21:03
好久了。这个问题好像还是没有解决。。。谢谢各位啦。

还没毕业么？

见楼下的回复。。。本楼重复了，请版主删掉？