Stata Code for Chapter Two evaluate.do
trueparam, delta(0) age(0) depress(0) predorig(5) h(0) numpat(1000) reseed(24) erseed(30)
display "age: "$age
display "depress: "$depress
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.ado.
(Stata ado program that performs the Monte Carlo simulations. This program is called by evaluate.do.) program trueparam, rclass
version 15.1
syntax, delta(numlist max=1 >=0) age(numlist integer >=0 <=1) depress(numlist integer >=0 <=1) predorig(integer) h(integer) ///
clear
[numpat(integer 1000) reseed(integer -1 ) erseed(integer -1)]
if `reseed'>=0&`erseed'>=0&`reseed'==`erseed'{
display as error "reseed has to be different from erseed."
exit, clear
}
*Seed for simulating random effects with drawnorm command if `reseed'<0 {
global reseed " "
}
else {
global reseed "seed(`reseed')"
}
*Seed for simulating model error terms
if `erseed'<0 {
global erseed " "
}
else {
global erseed "seed(`erseed')"
}
******************
global age=`age'
global depress=`depress'
global delta=`delta'
********************
global y=0.3365
6)) gives 0.3365
global PredOrig=`predorig'
// enter 1 if patient's age >65, 0 otherwise // enter 1 if patient had depression, 0 otherwisec
// The treatment target is <=6; The transformation log((6+1)/(11-
//PredOrig is the prediction origin (a time point).
global h=`h' //h=horizon
**********************
global NumPat=`numpat' //Enter number of simulated patients
* 4096 is the total number of possible combinations of true parameter values for a fixed value of delta.
* There are 12 model parameters that we need to calculate the true benefit and, therefore, 2^12=4096)
set matsize 4096
if $NumPat>4096 {
if $NumPat<=11000 set matsize $NumPat
else {
display as error "Number of simulated patients cannot be higher than 11000"
exit, clear
}
}
TrueParam_Module1
TrueParam_Module2
clear
quietly svmat Results, names(col)
save "Results_delta${delta}_age${age}_depress${depress}_PredOrig${PredOrig}_h${h}_NumPat${NumPat }.dta", replace
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) ********************
End
TrueParam_Module1.ado.
(Stata ado program used by trueparam.) set matsize 11000
**** This reads the estimates of the model reported in Table 1 of paper estimates use "Fitted_model.ster"
**** Covariance matrix
matrix D=e(cov_re)
global D11 D[1, 1] global D12 D[1, 2] global D13 D[1, 3] global D21 D[2, 1] global D22 D[2, 2]global D23 D[2, 3] global D31 D[3, 1] global D32 D[3, 2] global D33 D[3, 3]
**** Extract fixed effects matrix B=e(b)'
**** Fixed effects for pain score model global b4=_b[Marker:_M]
global b5=_b[Marker:_Mage_gt65] global b6=_b[Marker:_I_Mdepress_1]
global b7=_b[Marker:_Mtime]
global b8=_b[Marker:_I_MdX_Mtim_1] matrix b=( $b4 \ ///
$b5 \ /// $b6 \ /// $b7 \ /// $b8)
**** This gets the variance covariance matrix of fixed effects estimates
matrix VCe=e(V) // We will need the SE for fixed effects of the marker model in module 2. **** Variance of the error term for the pain score model in Table 1
global VarErr=e(var_eij)
**** Obtaining standard error of determinant of principal minor of D (eliminating 3rd row and 3rdcolumn) local D11 D[1, 1]
local D12 D[1, 2]
local D13 D[1, 3]
local D21 `D12'
local D22 D[2, 2]
local D23 D[2, 3]
local D31 `D13'
local D32 `D23'
local D33 D[3, 3]
**** Computation of principal minor of D and its SE
global pminor2=`D11'*`D22'-(`D21')^(2)
global SEpminor2=$pminor2*0.6 // We have to give a value since we cannot calculate the SE of principal minor
**** Computation of determinant of D and its standard error
Global detD=`D11'*(`D22'*`D33'-(`D32')^(2))-(`D21')^(2)*`D33'+2*`D21'*`D31'*`D32'- `D22'*(`D31')^(2)
global SEdetD=$detD*0.6
******************************************
End
TrueParam_Module2.ado.
(Stata ado program used by trueparam.) if $delta!=0 {
matrix Results=J(4096,2,0) local deltalist -$delta $delta
}
else {
matrix Results=J(1,2,0)
local deltalist 0 }
***********************
matrix colnames Results = RelBias Correlation 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' {
foreach delta7 of numlist `deltalist' {
foreach delta8 of numlist `deltalist' {
foreach delta9 of numlist `deltalist' {
foreach delta10 of numlist `deltalist' { foreach delta11 of numlist `deltalist' { foreach delta12 of numlist `deltalist' {
**** True fixed effects are computed
global b4True=$b4 + `delta1'*sqrt(VCe[4,4]) global b5True=$b5 + `delta2'*sqrt(VCe[5,5]) global b6True=$b6 + `delta3'*sqrt(VCe[6,6]) global b7True=$b7 + `delta4'*sqrt(VCe[7,7]) global b8True=$b8 + `delta5'*sqrt(VCe[8,8]) matrix bTrue=($b4True \ ///
$b5True \ /// $b6True \ /// $b7True \ ///
$b8True )
**** True variance covariance matrix of the random effects is computed
global D11True=$D11+`delta6'*0.228*0.2
global D12True=$D12+`delta7'*0.0931
global D13True=$D13+`delta8'*0.0942
global D21True=$D12True
local pminor2True=$pminor2+`delta9'*$SEpminor2 //The Variance-Covariance reparametrized to get a positive definite matrix
lobal D22True=(`pminor2True' +($D21True)^(2))/$D11True
global D23True=$D23+`delta10'*0.2008
global D31True=$D13True
global D32True=$D23True
matrix was
local detDTrue=$detD+`delta11'*$SEdetD //The Variance-Covariance matrix of random effects was reparametrized to get a positive definite matrix
global D33True=((`detDTrue'+$D11True*($D32True)^(2)-(2*$D21True*$D32True- $D22True*$D31True)*$D31True ) / ($D11True*$D22True-($D21True)^(2)))
matrix DTrue=( $D11True , $D12True , $D13True \ /// $D21True , $D22True , $D23True \ ///
$D31True , $D32True , $D33True ) *****************************************
**** True error variance is computed
global VarErrTrue=$VarErr+`delta12'* 0.38353*0.2
*TrueParam_Module3 simulates the random effects, creates design matrices of the appropriate size (according to PredOrig and h),
*and simulates the responses of patients, which are placed in matrix Y.
*The random effects plus their corresponding fixed effects are also saved in the database temporarily. display "Simulation `RowOfResults' for age=$age, depress=$depress,Delta=$delta, Prediction Origin=$PredOrig, Horizon=$h"
clear
TrueParam_Module3
*TrueParam_Module4 computes the BLUPs. They are saved in database.
TrueParam_Module4
*TrueParam_Module5 computes the empirical Bayesian predictors of benefits and true benefits (and save them temporarily in database)
TrueParam_Module5
matrix Results[`RowOfResults',1]=$RelBias
matrix Results[`RowOfResults',2]=$Correlation
local RowOfResults=`RowOfResults'+1
}
}
}
}
}
}
}
}
}
}
}
}
End
TrueParam_Module3.ado.
(Stata ado program used by trueparam.)
if $PredOrig==0|$PredOrig==1|$PredOrig==2|$PredOrig==3|$PredOrig==4 |$PredOrig==5 {
if $PredOrig==0 { matrix Z=(0, 1, 0 \ ///
1, 0, 0)
matrix X=(0, 0, 0, 1, 1*$age, 1*$depress, 0, 0*$depress \ ///
1, 1*$age, 1*$depress, 0, 0, 0, 0, 0)
local errors "eps01" }
if $PredOrig==1{ matrix Z=(0, 1, 0 \ ///
0, 1, 1 \ ///
1, 0, 0)
matrix X=(0, 0, 0, 1, 1*$age, 1*$depress, 0, 0*$depress \ ///
0, 0, 0, 1, 1*$age, 1*$depress, 1, 1*$depress \ ///
1, 1*$age, 1*$depress, 0, 0, 0, 0, 0) local errors "eps01 eps02"
}
if $PredOrig==2{ matrix Z=(0, 1, 0 \ ///
0, 1, 1 \ /// 0, 1, 2 \ /// 1, 0, 0)
matrix X=(0, 0, 0, 1, 1*$age, 1*$depress, 0, 0*$depress \ /// 0, 0, 0, 1, 1*$age, 1*$depress, 1, 1*$depress \ /// 0, 0, 0, 1, 1*$age, 1*$depress, 2, 2*$depress \ /// 1, 1*$age, 1*$depress, 0, 0, 0, 0, 0)
local errors "eps01 eps02 eps03" }
if $PredOrig==3{
matrix Z=(0, 1, 0 \ ///
0, 1, 1 \ /// 0, 1, 2 \ /// 0, 1, 3 \ /// 1, 0, 0)
matrix X=(0, 0, 0, 1, 1*$age, 1*$depress, 0, 0*$depress \ /// 0, 0, 0, 1, 1*$age, 1*$depress, 1, 1*$depress \ /// 0, 0, 0, 1, 1*$age, 1*$depress, 2, 2*$depress \ ///
0, 0, 0, 1, 1*$age, 1*$depress, 3, 3*$depress \ ///
1, 1*$age, 1*$depress, 0, 0, 0, 0, 0) local errors "eps01 eps02 eps03 eps04"
}
if $PredOrig==4{ matrix Z=(0, 1, 0 \ ///
0, 1, 1 \ /// 0, 1, 2 \ /// 0, 1, 3 \ /// 0, 1, 4 \ /// 1, 0, 0)
matrix X=(0, 0, 0, 1, 1*$age, 1*$depress, 0, 0*$depress \ ///
0, 0, 0, 1, 1*$age, 1*$depress, 1, 1*$depress \ ///
0, 0, 0, 1, 1*$age, 1*$depress, 2, 2*$depress \ /// 0, 0, 0, 1, 1*$age, 1*$depress, 3, 3*$depress \ /// 0, 0, 0, 1, 1*$age, 1*$depress, 4, 4*$depress \ ///
1, 1*$age, 1*$depress, 0, 0, 0, 0, 0) local errors "eps01 eps02 eps03 eps04 eps05"
}
if $PredOrig==5{ matrix Z=(0, 1, 0 \ ///
0, 1, 1 \ /// 0, 1, 2 \ /// 0, 1, 3 \ /// 0, 1, 4 \ /// 0, 1, 5 \ ///
1, 0, 0)
matrix X=(0, 0, 0, 1, 1*$age, 1*$depress, 0, 0*$depress \ ///
0, 0, 0, 1, 1*$age, 1*$depress, 1, 1*$depress \ /// 0, 0, 0, 1, 1*$age, 1*$depress, 2, 2*$depress \ /// 0, 0, 0, 1, 1*$age, 1*$depress, 3, 3*$depress \ /// 0, 0, 0, 1, 1*$age, 1*$depress, 4, 4*$depress \ ///
0, 0, 0, 1, 1*$age, 1*$depress, 5, 5*$depress \ /// 1, 1*$age, 1*$depress, 0, 0, 0, 0, 0)
local errors "eps01 eps02 eps03 eps04 eps05 eps06"
} }
else{
display as error "Values for predOrig should be 0,1,2,3,4, or 5"
exit, clear
}
matrix colnames Z=_D _M _Mtime
matrix colnames X=_D _Dage _Ddepress _M _Mage _Mdepress _Mtime _Mdepresstime *********************************************
if !(1<=$PredOrig+$h & $PredOrig+$h<=5) {
display as error "predorig+h should be in the interval [1,5]"
exit, clear
}
********************************************
*Variance covariance matrix of error terms
matrix RTrue=I(rowsof(Z)-1) * $VarErrTrue
*Simulation of pain scores Y
*First we simulate the patients
*The random effects are simulated.
*A particular value of vector (re_D re_M re_Mtime) correspond to one patient.
*The simulated random variables re_D re_M re_Mtime have mean zero
clear
quietly drawnorm re_D re_M re_Mtime , n($NumPat) cov(DTrue) $reseed
mkmat re_M re_Mtime, matrix(RanEff) // Each row of matrix RanEff corresponds to one simulated patient *Random coefficients are computed
*A random coefficient is what is usually called a random effect plus its corresponding fixed effect.
quietly generate re_M_c=$b4True+re_M
quietly generate re_Time_c=$b7True+re_Mtime
*The errors are simulated
quietly drawnorm `errors', n($NumPat) cov(RTrue) $erseed // errors are generated
mkmat `errors', matrix(Errors) // Each row of matrix Errors contains errors for corresponding patient in RanEff
*Columns of matrix MatbTrue will contain the true fixed effects repeatedly
matrix MatbTrue=bTrue
forvalues i=2/$NumPat {
matrix MatbTrue=MatbTrue,bTrue
}
*Matrix of responses is computed
*Each column of Y contains the simulated responses of the patient in corresponding column of RanEff' matrix Z2=Z[1..rowsof(Z)-1, 2..3]
matrix X2=X[1..rowsof(X)-1, 4..8]
matrix Y=Z2*RanEff'+X2*MatbTrue+Errors'
matrix colnames Y=Patient
matrix rownames Y=PainScore
end
2.7.6. TrueParam_Module4.ado.
(Stata ado program used by trueparam.)
*Variance covariance matrix of error terms
matrix R=I(rowsof(Z)) * $VarErr
matrix R[rowsof(Z),rowsof(Z)]=0
*Columns of matrix MatbGLS will contain the estimated fixed effects B repeatedly matrix MatbGLS=B
forvalues i=2/$NumPat {
matrix MatbGLS=MatbGLS,B
}
*Residual
matrix zero=0
matrix zero_pt=zero
forvalues i=2/$NumPat {
matrix zero_pt=zero_pt,zero
}
matrix Ynew=Y\zero_pt
matrix Res=Ynew-X*MatbGLS
matrix Res2=Res[1..rowsof(Res)-1, 1...]
matrix Res3=Res2\zero_pt
*The EB predictors are computed
matrix BLUP=D*Z'*inv(R+Z*D*Z')*Res3
matrix rownames BLUP=re1 re2 re3
*EB predictors are saved in database matrix BLUPT=BLUP'
svmat BLUPT, names(col)
end
TrueParam_Module5.ado.
(Stata ado program used by trueparam.)
*True benefit is computed
*(True variance of error is in global macro VarErrTrue)
quietly generate double TrueBenef=100*(normal( ($y-(re_M_c + $b5True*$age + $b6True*$depress +
re_Time_c *($PredOrig+$h) + $b8True*$depress*($PredOrig+$h) )
+ $b5True*$age + $b6True*$depress) ) /sqrt($VarErrTrue) )) label var TrueBenef "True benefit (x100)"
*Predicted benefit is computed
*(Computed with parameters in Table 1 of article)
) /sqrt($VarErrTrue) ) ///
- normal( ($y-(re_M_c
*(Variance of error from model in Table 1 is in global macro $VarErr)
*The EB predictors of the pain score model random effects are in re2 re3.
quietly generate double PredBenef=100*(normal( ($y-($b4+re2 + $b5*$age + $b6*$depress + ($b7+re3) *($PredOrig+$h) + $b8*$depress*($PredOrig+$h) ) ) /sqrt($VarErr) ) ///
($b4+re2 + $b5*$age + $b6*$depress) ) /sqrt($VarErr) )) label var PredBenef "Predicted benefit (x100)"
*Individual bias is computed
quietly generate double Bias=PredBenef-TrueBenef *Individual relative bias is computed
quietly summarize Bias, detail
- normal( ($y-
global MeanBias=r(mean) //Mean bias
quietly summarize TrueBenef, detail
global MeanTrueBenef=r(mean)
global RelBias=($MeanBias/$MeanTrueBenef)*100 // Relative bias quietly correlate TrueBenef PredBenef
global Correlation=r(rho)
end
留言
張貼留言