-
Notifications
You must be signed in to change notification settings - Fork 8
/
mrmvegger.ado
142 lines (115 loc) · 3.31 KB
/
mrmvegger.ado
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
*! version 0.1.0 23jun2020 Tom Palmer
program mrmvegger, eclass
version 9
local version : di "version " string(_caller()) ", missing :"
local replay = replay()
if replay() {
if _by() {
error 190
}
`version' Display `0', orientvar(`e(orientvar)') ///
n(`e(N)') np(`e(Np)') rmse(`e(phi)')
exit
}
syntax varlist(min=2) [aweight] [if] [in] [, ///
orient(integer 1) ///
Level(cilevel) ///
gxse(varlist numeric) ///
tdist ///
*]
if `orient' < 1 {
di as error "The orient() option cannot be negative"
exit 198
}
local callersversion = _caller()
// number of genotypes (i.e. rows of data used in estimation)
qui count `if' `in'
local k = r(N)
tokenize `"`varlist'"'
/*
varlist should be specified as:
1: gd beta
2: gp beta1
3: gp beta2
...
aw: =1/gdSE^2
*/
local outcome `1'
local npheno = wordcount("`varlist'") - 1
if `orient' > `npheno' {
di as error "The orient() option must be in the range of the number of phenotypes"
exit 198
}
local orientvar ``=`orient' + 1''
tempvar invvar gyse
qui gen double `invvar' `exp' `if' `in'
qui gen double `gyse' = 1/sqrt(`invvar') `if' `in'
tempvar gdtr eggercons
qui gen double `gdtr' = `1'*sign(``=`orient'+1'') / `gyse' `if'`in'
qui gen double `eggercons' = sqrt(`invvar') `if'`in'
local phenovarlist
local names
forvalues i = 1/`npheno' {
tempvar pheno`i'
qui gen double `pheno`i'' = ``=`i'+1'' * sign(``=`orient'+1'') / `gyse' `if'`in'
local phenovarlist "`phenovarlist' `pheno`i''"
local names `names' `outcome':``=`i'+1''
}
local names `names' `outcome':_cons
* fit model
tempname rmse
qui glm `gdtr' `phenovarlist' `eggercons' `if'`in', ///
nocons ///
level(`level') `options'
scalar `rmse' = sqrt(e(phi))
if `rmse' < 1 {
di as error "Residual standard error found to be:" as res scalar(`rmse'), ///
_n "This is less than 1.", ///
_n "Refitting model with residual variance constrained to 1."
local scale "scale(1)"
qui glm `gdtr' `phenovarlist' `eggercons' `if'`in', ///
nocons ///
level(`level') `options' `scale'
scalar `rmse' = sqrt(e(phi))
}
mat b = e(b)
mat V = e(V)
mat colnames b = `names'
mat rownames V = `names'
mat colnames V = `names'
ereturn post b V
ereturn local orientvar = "`orientvar'"
ereturn scalar N = `k'
ereturn scalar Np = `npheno'
ereturn scalar phi = scalar(`rmse')
if "`tdist'" != "" {
// use t-dist for ereturn display Wald test and CI limits
ereturn scalar df_r = `k' - `npheno' - 1
}
else {
// if df_r == . then Stata uses Normal dist
ereturn scalar df_r = .
}
* display estimates
Display , level(`level') orientvar(`orientvar') ///
n(`k') np(`npheno') rmse(`: di `rmse'')
ereturn local cmd "mrmvegger"
ereturn local cmdline `"mrmvegger `0'"'
end
program Display, rclass
syntax , [Level(cilevel)] orientvar(varname) N(integer) np(integer) rmse(real)
local orienttext : strlen local orientvar
local colstart = 79 - 31 - `orienttext'
di _n(1) _col(`colstart') as txt "MVMR-Egger model oriented wrt:", ///
as res "`orientvar'"
local nlength : strlen local n
local colstart = 79 - 22 - `nlength'
di _col(`colstart') as txt "Number of genotypes =", as res `n'
local nplength : strlen local np
local colstart = 79 - 23 - `nplength'
di _col(`colstart') as txt "Number of phenotypes =", as res `np'
di _col(47) as txt "Residual standard error =", as res %6.3f `rmse'
ereturn display, level(`level') noomitted
return add // r(table)
end
exit