I want to duplicate the resuls of sampleSelection package. However, my code does not work.
我想复制sampleSelection包的结果。然而,我的代码不起作用。
Data Set
library(mvtnorm)
set.seed(123)
n = 1000
sigmau <- 1; sigmav <- 1; rho <- 0.5
Sigma <- matrix(c(sigmau^2,rho*sigmau*sigmav,rho*sigmau*sigmav,sigmav^2),2,2)
x1 <- rnorm(n,0,1); x2 <- rnorm(n,0,1); z <- rnorm(n,0,1)
u = rmvnorm(n,c(0,0),Sigma)[,1]; v = rmvnorm(n,c(0,0),Sigma)[,2]
selection_latent <- 1 + x1 + z + v
selection <- as.numeric(selection_latent >= 0)
y <- ifelse(selection == 1, 1 + x1 + x2 + u, 0)
# Create the data frame 'df'
df_sel <- data.frame(id = 1:n,x1,x2,z,u = u,v,selection_latent,selection,y)
## Estimation
library(sampleSelection)
library(texreg)
sample_selection1 <- selection(selection = selection ~ 1 + x1 + z,outcome = y ~ 1 + x1 + x2, data = df)
screenreg(sample_selection1)
I wrote below code, but I receive the error message "Error in optim(par = parameters, fn = f_selection, selection = selection, :
function cannot be evaluated at initial parameters".
我编写了下面的代码,但收到错误消息“Error in OpTim(par=参数,fn=f_select,select=select,:函数不能在初始参数求值”)。
I would appriciate it if you find my error and modify the code. Thank you in advance.
如果你发现我的错误并修改代码,我会很感激。先谢谢你了。
## Sample Selection ##
Selection_Reg <- function(y_in, x_in, z_in, selection,df) {
x <- cbind(1,as.matrix(x_in))
z <- cbind(1,as.matrix(z_in))
lm1 <- lm(y~1+x1 +x2,data = df)
lm2 <- lm(y~1+x1 + z,data = df)
parameters <- c(1,1,as.numeric(lm1$coef),as.numeric(lm2$coef))
f_selection <- function(parameters,selection, y_in, x_in,z_in) {
sigma_u <- parameters[1]
rho <- parameters[2]
sigma_v <- 1
beta <- parameters[3:(ncol(x)+2)]
theta <- parameters[(ncol(x)+3):(length(parameters))]
xb <- x_in %*% beta
zt <- z_in %*% theta
z_adj <- (y_in-xb)/sigma_u
inside <- (zt + (rho*sigma_u/sigma_v)*(y_in-xb))/(((1-rho^2)*sigma_u^2)^(1/2))
log_lik <- (-log(sigma_u)+log(dnorm(z_adj))+log(pnorm(inside)))*selection +log(1 - pnorm(zt/sigma_v))*(1-selection)
return(-sum(log_lik))
}
result <- optim(par = parameters, fn = f_selection,selection = selection, y_in = y_in, x_in = x,z_in = z)
return(result)
}
Selection_Reg(y_in = df_sel$y,x_in = df_sel[,c("x1","x2")], z_in = df_sel[,c("x1","z")],selection = df_sel$selection,df = df_sel)
更多回答
优秀答案推荐
rho=1
in the initial values of the parameters passed to optim()
is the main issue, since it makes the variable inside
undefined. Changing rho=0.5
in the function Selection_Reg()
(same initialization as in the first case) makes the code work:
传递给OpTim()的参数的初始值中的RHO=1是主要问题,因为它使内部变量未定义。在函数SLECTION_REG()中更改RHO=0.5(与第一种情况相同的初始化)可以使代码正常工作:
parameters <- c(1,0.5,as.numeric(lm1$coef), as.numeric(lm2$coef))
Now, running optim()
converges to the following (MLE) values:
现在,运行OpTim()收敛到以下(MLE)值:
Selection_Reg(y_in = df_sel$y,
x_in = df_sel[,c("x1","x2")],
z_in = df_sel[,c("x1","z")],
selection = df_sel$selection, df = df_sel)
# $par
# [1] 0.9460186 0.2447177 0.8799579 1.0583746 0.9466844 0.9414059 0.9505922 0.8971660
# $value
# [1] 1351.807
# $counts
# function gradient
# 501 NA
# $convergence
# [1] 1
Where in the first case the output is the following:
其中,在第一种情况下,输出如下:
screenreg(sample_selection1)
#
#============================
# maximum
#----------------------------
#S: (Intercept) 0.93 ***
# (0.06)
#S: x1 1.04 ***
# (0.07)
#S: z 0.98 ***
# (0.07)
#O: (Intercept) 0.99 ***
# (0.07)
#O: x1 0.99 ***
# (0.05)
#O: x2 0.96 ***
# (0.04)
#sigma 0.99 ***
# (0.03)
#rho -0.13
# (0.16)
#----------------------------
#AIC 2718.13
#BIC 2757.40
#Log Likelihood -1351.07
#Num. obs. 1000
#Censored 293
#Observed 707
#============================
#*** p < 0.001; ** p < 0.01; * p < 0.05
Now, compare MLE values of sigma, rho, beta, theta parameters, obtained with the above two approaches.
现在,比较上述两种方法得到的sigma、Rho、beta、theta参数的最大似然估计值。
更多回答
Dear Dey, Thanks a lot! I really appriciate your help :)
亲爱的亲爱的,非常感谢!我真的很感谢你的帮助:)
我是一名优秀的程序员,十分优秀!