将离散选择实验样本量计算的R代码转换为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)
matrix parameters = parameters'  // 转置为列向量，适配矩阵乘法

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

// 加载设计矩阵数据
import delimited "~/Downloads/Design_matrix_Illustration_DCE_Osteoporosis_treatment.txt", clear
mkmat v1 v2 v3 v4 v5 v6 v7 v8, matrix(design)

// 初始化信息矩阵
matrix info_mat = J(ncoefficients, ncoefficients, 0) 

// 计算exp(设计矩阵 × 参数向量)：修正点积运算
mata:
    design = st_matrix("design")
    params = st_matrix("parameters")
    exputilities = exp(design * params)
    st_matrix("exputilities", exputilities)
end

// 循环处理每个选择集，累积信息矩阵
forvalues k_set = 1/`nchoices' {
    local start = (`k_set'-1)*`nalts' + 1
    local end = `k_set'*`nalts'
    
    // 提取当前选择集的设计矩阵行和效用值
    matrix design_sub = design[`start'..`end', .]
    matrix exp_sub = exputilities[`start'..`end', .]
    
    // 计算选择概率p_set
    scalar sum_exp = rowtotal(exp_sub)'  // 转置为标量
    matrix p_set = exp_sub / sum_exp
    
    // 构造对角矩阵P和外积pp'
    mata:
        p = st_matrix("p_set")
        p_diag = diag(p)
        pp_t = p * p'
        middle_term = p_diag - pp_t
        st_matrix("middle_term", middle_term)
    end
    
    // 计算当前选择集对信息矩阵的贡献
    matrix full_term = design_sub' * middle_term * design_sub
    matrix info_mat = info_mat + full_term
}

// 计算方差-协方差矩阵（信息矩阵的逆）
matrix sigma_beta = inv(info_mat)

// 计算每个系数所需样本量
matrix effectsize = parameters  // 使用参数值作为效应量
mata:
    z_beta = st_scalar("z_one_minus_beta")
    z_alpha = st_scalar("z_one_minus_alpha")
    sigma_diag = diag(st_matrix("sigma_beta"))
    effect = abs(st_matrix("effectsize"))
    N = ((z_beta + z_alpha) :* sqrt(sigma_diag) :/ effect) :^ 2
    st_matrix("N", N)
end

// 展示结果
matrix list N, noheader

关键步骤解释 #

1. 修正点积运算 #

原R代码中design%*%parameters是矩阵与向量的点积，Stata需将参数向量转置为列向量，再通过Mata完成矩阵乘法（Mata的矩阵运算逻辑更贴近R）。

2. 选择集循环处理 #

• 用forvalues遍历每个选择集，通过索引范围提取对应行的设计矩阵和效用值
• 计算选择概率p_set：当前选择集内各选项效用值的占比
• 构造中间项P - pp'：通过Mata的diag()生成对角矩阵，再计算外积并做差

3. 信息矩阵累积 #

每个选择集的贡献项full_term通过矩阵转置和乘法计算后，累加到初始的零矩阵info_mat中。

4. 样本量计算 #

• 对信息矩阵求逆得到方差-协方差矩阵sigma_beta
• 提取对角线元素（各系数的方差），代入样本量公式计算每个系数所需的最小样本量

内容的提问来源于stack exchange，提问作者Obvious Nchimunya Chilyabanyam