|
1 | 1 | """
|
2 | 2 | lrt(model1::VCModel, model0::VCModel)
|
3 | 3 |
|
4 |
| -Perform a variation of the likelihood ratio test for univariate variance components models as in |
5 |
| -Molenberghs and Verbeke 2007 with model1 and model0 being the full and nested models, respectively. |
| 4 | +Perform a variation of the likelihood ratio test for univariate variance components models |
| 5 | +as in Molenberghs and Verbeke 2007 with `model1` and `model0` being the full and nested |
| 6 | +models, respectively. |
6 | 7 | """
|
7 | 8 | function lrt(
|
8 | 9 | model1 :: VCModel,
|
9 | 10 | model0 :: VCModel
|
10 | 11 | )
|
11 | 12 | df = length(model1.V) - length(model0.V)
|
12 |
| - @assert df > 0 |
13 |
| - @assert size(model0.Σ[1], 1) == 1 |
14 |
| - @assert size(model1.Σ[1], 1) == 1 |
| 13 | + @assert df > 0 "models should be nested" |
| 14 | + @assert size(model1.Σ[1], 1) == 1 && size(model0.Σ[1], 1) == 1 "models should be univariate" |
15 | 15 | lrt = 2 * (model1.logl[1] - model0.logl[1])
|
16 | 16 | coefs = [2.0^-df * binomial(df, i) for i in 1:df]
|
17 | 17 | ps = [ccdf(Chisq(i), lrt) for i in 1:df]
|
18 |
| - sum(coefs .* ps) |
| 18 | + dot(coefs, ps) |
19 | 19 | end
|
20 | 20 |
|
21 | 21 | """
|
22 | 22 | h2(model::VCModel)
|
23 | 23 |
|
24 |
| -Calculate heritability estimates and their standard errors, assuming that all variance components |
25 |
| -capture genetic effects except the last term. Also return total heritability from sum of individual |
26 |
| -contributions and its standard error. |
| 24 | +Calculate heritability estimates and their standard errors, assuming that all variance |
| 25 | +components capture genetic effects except the last term. Also return total heritability from |
| 26 | +sum of individual contributions and its standard error. |
27 | 27 | """
|
28 | 28 | function h2(model::VCModel)
|
29 | 29 | m, d = length(model.Σ), size(model.Σ[1], 1)
|
@@ -52,16 +52,6 @@ function h2(model::VCModel)
|
52 | 52 | h2s, ses
|
53 | 53 | end
|
54 | 54 |
|
55 |
| -function findvar(d::Int) |
56 |
| - s, r = ◺(d), d |
57 |
| - idx = ones(Int, d) |
58 |
| - for j in 2:length(idx) |
59 |
| - idx[j] = idx[j - 1] + r |
60 |
| - r -= 1 |
61 |
| - end |
62 |
| - idx |
63 |
| -end |
64 |
| - |
65 | 55 | """
|
66 | 56 | rg(model::VCModel)
|
67 | 57 |
|
@@ -94,3 +84,13 @@ function rg(model::VCModel)
|
94 | 84 | [copytri!(ses[i], 'L') for i in 1:m]
|
95 | 85 | rgs, ses
|
96 | 86 | end
|
| 87 | + |
| 88 | +function findvar(d::Int) |
| 89 | + s, r = ◺(d), d |
| 90 | + idx = ones(Int, d) |
| 91 | + for j in 2:length(idx) |
| 92 | + idx[j] = idx[j - 1] + r |
| 93 | + r -= 1 |
| 94 | + end |
| 95 | + idx |
| 96 | +end |
0 commit comments