Stata Code for Chapter One: evaluate_disorganized.do

 trueparam_disorganized, delta(0.5) canna01(0) predorig(3) h(1) numpat(1000) reseed(24)

display "canna01: "$canna01

display "Prediction origin: t="$PredOrig

display "Prediction horizon: h="$h

display "delta= " $delta

display "Median of relative biases: " r(MedianRelBias)

display "Minimum of relative biases: " r(MinRelBias)

display "Maximum of relative biases: " r(MaxRelBias)

display "Median of correlations between predicted and true transformed benefits: " r(MedianCorr) display "Minimum of correlations: " r(MinCorr)

display "Maximum of correlations: " r(MaxCorr)

trueparam_disorganized.ado.

(Stata ado program that performs the Monte Carlo simulations. This program is called by evaluate_disorganized.do)

program trueparam_disorganized, rclass

version 15.1

syntax, delta(numlist max=1 >=0) canna01(numlist integer >=0 <=1) predorig(integer) h(integer) ///

[numpat(integer 1000) reseed(integer -1 )]

clear

*Seed for simulating random effects with drawnorm command if `reseed'<0 {

global reseed " "


}

else {

global reseed "seed(`reseed')"

}

******************

global canna01=`canna01' // enter 1 if patient used cannabis

*******************

global PredOrig=`predorig' //PredOrig is the prediction origin (a time point).

global h=`h' //h=horizon; enter a negative number or 0 for retrospective measurement of benefits; **********************

global NumPat=`numpat' //Enter number of simulated patients

*******************

global delta=`delta' // Enter distance from a true parameter to parameter estimate in Table 1 in standard error units

**********************

* 64 is the total number of possible combinations of true parameter values for a fixed value of delta.

* There are 6 model parameters and, therefore, 2^ 6 = 64

set matsize 64

if $NumPat>64 {

if $NumPat<=11000 set matsize $NumPat

else {

display as error "Number of simulated patients cannot be higher than 11000"

exit, clear

}

}

TrueParam_Module1_disorganized


TrueParam_Module2_disorganized

clear

quietly svmat Results, names(col) // convert the matrix into a data set

save "Results_delta${delta}_canna01${canna01}_PredOrig${PredOrig}_h${h}_NumPat${NumPat}.dta", replace

quietly summarize MeanBias,detail

return scalar MedianMeanBias=r(p50)

return scalar MinMeanBias=r(min)

return scalar MaxMeanBias=r(max)

quietly summarize SDBias, detail

return scalar MedianSDBias=r(p50)

return scalar MinSDBias=r(min)

return scalar MaxSDBias=r(max)

quietly summarize RelBias, detail

return scalar MedianRelBias=r(p50)

return scalar MinRelBias=r(min)

return scalar MaxRelBias=r(max)

quietly summarize Correlation, detail

return scalar MedianCorr=r(p50)

return scalar MinCorr=r(min)

return scalar MaxCorr=r(max)

quietly summarize MeanTrueBenef, detail

return scalar MedianMeanTrueBenef=r(p50)

return scalar MinMeanTrueBenef=r(min)

return scalar MaxMeanTrueBenef=r(max)


********************

end

TrueParam_Module1_disorganized.ado.

(Stata ado program used by trueparam_disorganized.) set more off

set matsize 11000

*This reads the estimates of the model reported in Table 1 of paper estimates use "Fitted_model"

*This gets the variance covariance matrix of estimates

matrix VCe=e(V) ********************************************

*Vector of estimated fixed effects from Model in Table 1 is created global b1=_b[dis_lt4:_cons]

global b2=_b[dis_lt4:canna01]

global b3=_b[dis_lt4:pt1]

global b4=_b[dis_lt4:pt2]

global b5=_b[dis_lt4:pt3]

matrix bGLS=($b1 \ ///

$b2 \ /// $b3 \ ///

$b4 \ ///

$b5)

*The variance-covariance matrix D of random effects from model in Table 1 is created global D11=_b[/var(_cons[id])]

matrix D=$D11

end


TrueParam_Module2_disorganized.ado.

(Stata ado program used by trueparam_disorganized.) if $delta!=0 {

matrix Results=J(64,5,0) //2^6=64 local deltalist -$delta $delta

}

else {

matrix Results=J(1,5,0) local deltalist 0

}

***********************

matrix colnames Results = MeanBias SDBias RelBias Correlation MeanTrueBenef local RowOfResults=1

foreach delta1 of numlist `deltalist' {

foreach delta2 of numlist `deltalist' {

foreach delta3 of numlist `deltalist' {

foreach delta4 of numlist `deltalist' {

foreach delta5 of numlist `deltalist' {

foreach delta6 of numlist `deltalist' {

*True fixed effects are computed

global b1True=$b1 + `delta1'*sqrt(VCe[5,5])

global b2True=$b2 + `delta2'*sqrt(VCe[1,1])

global b3True=$b3 + `delta3'*sqrt(VCe[2,2])

global b4True=$b4 + `delta4'*sqrt(VCe[3,3])

global b5True=$b5 + `delta5'*sqrt(VCe[4,4])

matrix bTrue=($b1True \ ///


$b2True \ /// $b3True \ /// $b4True \ /// $b5True)

matrix list bTrue

*True variance covariance matrix is computed

global D11True=$D11+`delta6'*sqrt(VCe[6,6])

matrix DTrue=$D11True

*****************************************

display "Simulation `RowOfResults' for canna01=$canna01, Delta=$delta, Prediction Origin=$PredOrig, Horizon=$h"

clear

TrueParam_Module3_disorganized

matrix Results[`RowOfResults',1]=$MeanBias

matrix Results[`RowOfResults',2]=$SDBias

matrix Results[`RowOfResults',3]=$RelBias

matrix Results[`RowOfResults',4]=$Correlation

matrix Results[`RowOfResults',5]=$MeanTrueBenef

local RowOfResults=`RowOfResults'+1

}

}

}

}

}

}

End

/* variance of intercept */


TrueParam_Module3_disorganized.ado.

(Stata ado program used by trueparam_disorganized.) *Design matrix Z for random effects

matrix A=J(7, 1, 1)

matrix pt1=J(7, 1, 0)

forvalues i=1/7{

matrix pt1[`i', 1]=P[1,1]*(`i'-1) + P[1,2]*(`i'-1)^2 + P[1,3]*(`i'-1)^3 + P[1,4]

}

matrix Z=A

matrix colnames Z=Intercept ********************************************** *Design matrix X for fixed effects

matrix A=J(7, 1, 1)

matrix B=J(7, 1, $canna01)

matrix pt2=J(7, 1, 0)

forvalues i=1/7{

matrix pt2[`i', 1]=P[2,1]*(`i'-1) + P[2,2]*(`i'-1)^2 + P[2,3]*(`i'-1)^3 + P[2,4] }

matrix pt3=J(7, 1, 0) forvalues i=1/7{

matrix pt3[`i', 1]=P[3,1]*(`i'-1) + P[3,2]*(`i'-1)^2 + P[3,3]*(`i'-1)^3 + P[3,4] }

matrix X=A,B,pt1,pt2,pt3

matrix colnames X=Intercept canna01 ptime1 ptime2 ptime3 **********************************************


if $PredOrig==0|$PredOrig==1|$PredOrig==2|$PredOrig==3|$PredOrig==4|$PredOrig==5 |$PredOrig==6 {

if $PredOrig==0 { matrix Z=Z[1..1,1...] // matrix X=X[1..1,1...]

}

if $PredOrig==1{ matrix Z=Z[1..2,1...] matrix X=X[1..2,1...] }

if $PredOrig==2{ matrix Z=Z[1..3,1...] matrix X=X[1..3,1...] }

if $PredOrig==3{ matrix Z=Z[1..4,1...] matrix X=X[1..4,1...] }

if $PredOrig==4{ matrix Z=Z[1..5,1...] matrix X=X[1..5,1...]

}

if $PredOrig==5{

matrix Z=Z[1..6,1...] matrix X=X[1..6,1...]


}

if $PredOrig==6{ matrix Z=Z[1..7,1...] matrix X=X[1..7,1...] }

}

else{

display as error "Values for predOrig should be 0,1,2,3,4,5,or 6"

exit, clear

}

*********************************************

if !(1<=$PredOrig+$h & $PredOrig+$h<=6) {

display as error "predorig+h should be in the interval [1,6]"

exit, clear

}

********************************************

*First we simulate the patients

*The simulated random intercept LambdaR has mean zero

clear

quietly drawnorm LambdaR, n($NumPat) cov(DTrue) $reseed // DTrue was created in Module2 mkmat LambdaR, matrix(RanEff) // Each row of matrix RanEff corresponds to a set of random effects (for intercept and ptime1) for one simulated patient

*Columns of matrix MatbTrue will contain the true fixed effects repeatedly

matrix MatbTrue=bTrue

forvalues i=2/$NumPat {


matrix MatbTrue=MatbTrue,bTrue

// each column has a vector of bTrue

}

*Calculate linear predictor

matrix RanEfft=RanEff'

matrix ZR=Z*RanEfft

matrix XB=X*MatbTrue

matrix ata=ZR+XB

*For each element of anta, calculate p. This will be the predicted P, which can be used to calculate the predicted benefit as well as to simulate y.

matrix p=J(rowsof(ata), colsof(ata), 0) // number of time points by number of patients

matrix y=J(rowsof(ata), colsof(ata), 0)

/* begin loop */

local i=1 // i for time

while `i'<=rowsof(ata){

local j=1 // j for patients while `j'<=colsof(ata){

matrix p[`i',`j']=exp(ata[`i',`j'])/(1+exp(ata[`i',`j']))

matrix y[`i',`j']=rbinomial(1, p[`i',`j'])

local j = `j' + 1 }

local i = `i' + 1 }

/* end loop */

*Create dataset xy with y and Xs for 1000 patients in rows by appending matrices. Include patient id


matrix X1=X

forvalues i=2/$NumPat {

matrix X1=X1\X }

matrix y1=y[1..., 1] forvalues j=2/$NumPat {

matrix y2=y[1..., `j']

matrix y1=y1\y2 }

matrix id=J(rowsof(ata), 1, 1) forvalues j=2/$NumPat {

matrix id2=J(rowsof(ata), 1, `j')

matrix id=id\id2 }

matrix xy=y1, X1, id

matrix colnames xy=dis_lt4 Intercept canna01 pt1 pt2 pt3 id

// convert to dataset

clear

svmat xy, names(col)

predict PrRanEff* , reffects

duplicates drop id, force

keep PrRanEff* id

svmat RanEff, names(col)

spearman PrRanEff1 LambdaR ******************************************************************************** *Compute predicted benefit and true benefit for each patient


quietly generate ppred2=1/(1+ exp(-($b1 + PrRanEff1 + $b2*$canna01 + $b3 * pt1[$PredOrig+$h+1, 1] + /// $b4 * pt2[$PredOrig+$h+1, 1] + $b5 * pt3[$PredOrig+$h+1, 1] ))) quietly generate ppred1=1/(1+ exp(-($b1 + PrRanEff1 + $b2*$canna01 + $b3* pt1[1, 1] + ///

$b4 * pt2[1, 1] + $b5 * pt3[1, 1] )))

quietly generate PredBenef=ppred2-ppred1

quietly generate ptrue2=1/(1+ exp(-($b1True + LambdaR + $b2True*$canna01 + $b3True * pt1[$PredOrig+$h+1, 1] + ///

$b4True * pt2[$PredOrig+$h+1, 1] + $b5True *

pt3[$PredOrig+$h+1, 1] )))

quietly generate ptrue1=1/(1+ exp(-($b1True + LambdaR + $b2True*$canna01 + $b3True * pt1[1, 1]+ ///

$b4True * pt2[1, 1] + $b5True * pt3[1, 1] )))

quietly generate TrueBenef=ptrue2-ptrue1

keep PredBenef TrueBenef ******************************************************************************** *Individual bias is computed

quietly generate double Bias=PredBenef-TrueBenef

*Individual relative bias is computed

quietly summarize Bias, detail

global MeanBias=r(mean) //Mean bias

global SDBias=sqrt(r(Var)) // SD of bias

quietly summarize TrueBenef, detail

global MeanTrueBenef=r(mean)

global RelBias=($MeanBias/$MeanTrueBenef)*100 // Relative bias ******************************************************************************** *Computation of correlation between predicted and true benefit


quietly spearman TrueBenef PredBenef global Correlation=r(rho)

end

留言

張貼留言

這個網誌中的熱門文章

可轉移性、普遍性、代表性和外部有效性

 ### 可轉移性（Transportability） - **定義：** 可轉移性指的是將一個研究中得到的因果效應估計從一個人群（研究人群）應用到另一個不同的人群（目標人群）的能力。它涉及特定的統計方法來調整人群之間的差異，並特別關注為了合理和實施此類應用所需的假設和方法。 - **用途：** 當研究人員想知道在一個設置或人群中觀察到的干預效果是否在不同的設置或人群中保持時使用。 ### 普遍性（Generalizability） - **定義：** 普遍性是指研究發現可以應用到更廣泛人群、設置、時間或未直接研究的測量指標的程度。它是一個比可轉移性更廣泛的概念，更多地關注研究結果在特定研究條件之外的適用性。 - **用途：** 在評估研究結果是否適用於明確研究之外的更廣泛背景時考慮，而不必然採用特定的統計方法來進行此類應用。 ### 代表性（Representativeness） - **定義：** 代表性指的是研究樣本如何緊密地反映其所來自的總體或目標人群的特徵。它涉及抽樣過程和樣本在多大程度上體現了更廣泛人群的特徵。 - **用途：** 這一概念在研究設計階段至關重要，影響研究結果能夠被泛化到總體的程度。 ### 外部有效性（External Validity） - **定義：** 外部有效性是指研究結果可以被泛化或應用到原始研究設置、人群、治療變量和結果之外的背景的程度。它包括普遍性和可轉移性的概念，但範圍更廣，包括隨時間和在不同條件下應用研究發現的適用性。 - **用途：** 它涉及到研究發現超出特定研究範圍的整體適用性和相關性，考慮諸如時間、設置和人群等因素。 **關鍵差異和關係：** - **可轉移性與普遍性：** 可轉移性是普遍性的一個子集，特別關注使用統計調整將結果應用到不同人群的能力。普遍性是一個更廣泛的概念，不必然涉及此類調整。 - **代表性與外部有效性：** 代表性樣本可以增強研究的外部有效性，但外部有效性也取決於其他因素，如研究設計和執行。 - **外部有效性包含：** 普遍性和可轉移性，以及其他考慮因素，以便將研究發現應用於原始背景之外。 總之，雖然所有這些概念都涉及到將研究發現應用於原始研究條件之外，但它們每一個都關注此應用的不同方面，從用於調整不同人群之間差異的統計方法（可轉移性）到研究發現的整體適用性（外部有效性）。 - ...
閱讀完整內容

頻率學派 vs 貝氏學派

機率是重複事件長期發生的頻率 vs 主觀上某事件發生的不確定性（比較像是一般人的觀念）。 傳統 vs 新 簡單vs 複雜 人工計算可行 vs 不可行 母數（參數）是未知的常數（沒有機率可言） vs 未知的隨機變數（有機率） 樣本是隨機的變數 vs 已經發生的固定數 客觀（只依據本次試驗的結果）vs 主觀（依據先驗信念及本次試驗的結果而更新成後驗信念）。 H: 假設，D: 資料。H0: 虛無假設（參數或母數 θ 等於 0），H1: 對立假設，p(H1)=1-p(H0)。 一個棒球比賽的觀眾關心的是「機率」p(D|H)，亦即預測「這個左打者擊中球的機率是多少？」；當他看見一名打者擊中了球以後，計算這一名打者是一個左打者的機率是多少則是屬於「逆機率」p(H|D)。 醫生關心的是「逆機率」p(H|D)：當看見某症候、檢驗、檢查 D時，用貝氏定理計算病人有某疾病 H 的機率。 估計 p 值（H0 為真時，看見 D 或比它更極端的機率 p(D|H0)) vs 估計後驗機率 p(H1|D)：p(D|H) ≠ p(H|D)。 95% 信賴區間包含母數的機率是 0% 或 100% vs 95% 置信區間包含母數的機率是 95%。 頻率學派：當 p < α（例如 0.05）時拒絕 H0，當 p ≥ α 時無法拒絕 H0。但是我們無法知道 p(H0|D)、p(H1|D)，而且無論 p 值多大都無法接受 H0。 貝氏定理： 1. H1 的後驗勝算 p(H1|D)/(p(H0|D))= 先驗勝算 x 似然比或貝氏因子 p(D|H1)/p(D|H0)，勝算=機率/(1-機率），機率=勝算/(勝算+1)。 2. 後驗機率 p(H|D）= 先驗機率 p(H) x p(D/H)/p(D) p(D): p(D/H0) + p(D/H1) RCT 時，頻率學派的 α 是 H0 為真（無效）時，你卻斷言為假（有效）的機率：p(斷言有效|無效)。傳統上它被稱為「第一型錯誤率」，其實它應該是「第一型斷言率」，而不是治療無效你卻以為有效的錯誤機率：貝氏學派的 1 - p(有效|資料，先驗機率）。後者才是我們所關心的。 ESP and the significance of significance 22 https://understandinguncertainty.org/node/1286 爲什麼我...
閱讀完整內容

貝氏分析計算器

貝氏分析計算器 ：20250130 by 顧進裕, Claude-3.5-Sonnet 平均值與標準差計算器 貝氏定理：後驗分布 p(H|D) = 先驗分布 p(H) × 似然分布 p(D|H) / 邊際分布 p(D) 似然分布：觀察到的平均值與標準誤 se 後驗分布：95% 可信區間（credible interval） 1. 等尾可信區間：適用於對稱性分布，可以跟 95% CI 比較 2. 最高後驗密度：適用於偏斜或多峰分布 3. 分位數可信區間：適用於單峰的對稱性分布 後驗預測分布：95% 預測區間，p 值靠近 0.5 時表示模型是可靠的 假設： • log(OR), log(RR), log(HR) 都是常態分布 • 先驗 log(OR) = log(RR) = log(HR) = 0（即先驗 OR = RR = HR = 1） • log(先驗反應比率)=log(mean 1/mean2) • 大部分的 OR, RR, HR, log(先驗反應比率) 介於 0.5～2 之內  • 先驗標準平均差 SMD =平均差/標準差 sd= 0 • 大部分的 SMD 介於 -1～1（亦即二組相差一個標準差）之內 多重先驗敏感性分析： • 強先驗: gaussian(0, 0.2)，95% confidence interval（信賴區間） -0.4, 0.4, 亦即 OR, RR, HR, 先驗反應比率 0.7, 1.5 • 內定先驗: gaussian(0, 0.42)，95% CI（信賴區間） -0.84, 0.84, 亦即 OR, RR, HR, 先驗反應比率 0.4, 2.7 • 弱先驗: gaussian(0, 1)，95% CI（信賴區間） -4, 4, 亦即 OR, RR, HR, 先驗反應比率 0.14, 7.4
閱讀完整內容