将 R 代码翻译成 Stata 代码以计算离散选择实验的样本量

Question

我正在尝试计算离散选择实验的样本量，但我发现的文献和代码是用 R 编写的。但是我正在尝试在 Stata 中复制代码。

test_alpha=0.05
z_one_minus_alpha<-qnorm(1-test_alpha)

test_beta=0.2
z_one_minus_beta<-qnorm(1-test_beta)

parameters<-c(1.23 , -0.31 , -0.21 , -0.44 , 0.028 , -1.10 , -0.04 , -0.0015)

ncoefficients=8
nalts=3
nchoices=16

# load the design information
design<-as.matrix(read.table("~/Downloads/Design_matrix_Illustration_DCE_Osteoporosis_treatment.txt", header=FALSE));

#compute the information matrix
# initialize a matrix of size ncoefficients by ncoefficients filled with zeros.
info_mat=matrix(rep(0,ncoefficients*ncoefficients), ncoefficients, ncoefficients) 
# compute exp(design matrix times initial parameter values) 
exputilities=exp(design%*%parameters)
# loop over all choice sets
for (k_set in 1:nchoices) {
  # select alternatives in the choice set
  alternatives=((k_set-1)*nalts+1) : (k_set*nalts)
  # obtain vector of choice shares within the choice set
  p_set=exputilities[alternatives]/sum(exputilities[alternatives])
  # also put these probabilities on the diagonal of a matrix that only contains zeros
  p_diag=diag(p_set)
  # compute middle term P-pp’
  middle_term<-p_diag-p_set%o%p_set
  # pre- and postmultiply with the Xs from the design matrix for the alternatives in this choice set
  full_term<-t(design[alternatives,])%*%middle_term%*%design[alternatives,]
  # Add contribution of this choice set to the information matrix
  info_mat<-info_mat+full_term 
} # end of loop over choice sets
#get the inverse of the information matrix (i.e., gets the variance-covariance matrix)
sigma_beta<-solve(info_mat,diag(ncoefficients)) 

# Use the parameter values as effect size. Other values can be used here.
effectsize<-parameters
# formula for sample size calculation is n>[(z_(beta)+z_(1-alpNha))*sqrt(S??)/delta]^2 
N<-((z_one_minus_beta + z_one_minus_alpha)*sqrt(diag(sigma_beta))/abs(effectsize))^2
# Display results (required sample size for each coefficient)
N

我在 Stata 中尝试了以下代码，但是在我需要 Stata 中的点积并运行最后一个循环的地方我得到了堆栈：

scalar test_alpha = 0.05
scalar z_one_minus_alpha = invnormal(1-test_alpha)

scalar test_beta = 0.2
scalar z_one_minus_beta = invnormal(1-test_beta)

matrix parameters = (1.23, -0.31, -0.21, -0.44, 0.028, -1.10, -0.04, -0.0015)

scalar ncoefficients = 8
scalar nalts = 3
scalar nchoices = 16

// load the design information
import delimited "~/Downloads/Design_matrix_Illustration_DCE_Osteoporosis_treatment.txt", clear

/*adding a column of 1s for the coeffficients*/
*gen v0 =1
*order v0 v1 v2 v3 v4 v5 v6 v7 v8

mkmat v1 v2 v3 v4 v5 v6 v7 v8,matrix(design)

//compute the information matrix
// initialize a matrix of size ncoefficients by ncoefficients filled with zeros.
matrix info_mat = J(ncoefficients, ncoefficients, 0) 

// compute exp(design matrix times initial parameter values) 

matrix A = parameters#design

mata : st_matrix("exputilities", exp(st_matrix("A")))

mat li exputilities

到目前为止，我已经编写了上面的代码。