/
admixture-weight-estimation.Rmd
167 lines (140 loc) · 9.17 KB
/
admixture-weight-estimation.Rmd
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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
---
title: "Estimation of unknown elements in admixture models"
author: "Xavier Milhaud"
bibliography: references.bib
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Estimation of unknown elements in admixture models}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
warning = FALSE,
message = FALSE
)
```
```{r setup}
library(admix)
```
We remind that a random variable $X$ following an admixture distribution has cumulative distribution function (cdf) $L$ given by
$$L(x) = pF(x) + (1-p)G(x), \qquad x \in \mathbb{R},$$
where $G$ is a mixture component whose distribution is perfectly known, whereas $p$ and $F$ are unknown.\\
In this setting, if no parametric assumption is made on the unknown component distribution $F$, the mixture is considered as a semiparametric mixture. For an overview on semiparametric extensions of finite mixture models, see [@Xiang19].
# Estimation of the unknown component weight in an admixture model
The mixture weight $p$ of the unknown component distribution can be estimated using diverse techniques depending on the assumptions made on the unknown cdf $F$, among which the ones discussed in the sequel:
* the estimator provided by Bordes and Vandekerkhove, see [@BV10];
* the estimator provided by Patra and Sen, [@Patra];
* the estimator provided by Inversion - Best Matching method, [@MilhaudIBM21].
All the estimation methods presented hereafter can be performed using one single generic function for estimation with appropriate arguments, the so-called $admix_estim$ function.
## The one-sample case
Many works studied the estimation of the unknown proportion in two-component admixture models. Among them, seminal papers are [@BDV06] and [@BMM06]. These papers are closely connected to the paper by [@BV10], where an asymptotic normal estimator is provided for the unknown component weight.
### Case of symmetric unknown density
In this case, we use the Bordes and Vandekerkhove estimator, see [@BV10].
```{r}
## Simulate data:
list.comp <- list(f = 'norm', g = 'norm')
list.param <- list(f = list(mean = 3, sd = 0.5), g = list(mean = 0, sd = 1))
data1 <- rsimmix(n = 800, unknownComp_weight = 0.7, list.comp, list.param)[['mixt.data']]
## Perform the estimation of parameters in real-life:
list.comp <- list(f = NULL, g = 'norm')
list.param <- list(f = NULL, g = list(mean = 0, sd = 1))
BVdk_estimParam(data1, method = 'L-BFGS-B', list.comp, list.param)
```
Here we look at the first listed parameter, which corresponds to the estimate of the unknown component weight.
### Other cases
In this case, we use the Patra and Sen estimator, see [@Patra].
```{r}
## Simulate data:
list.comp <- list(f = 'norm', g = 'norm')
list.param <- list(f = list(mean = 3, sd = 0.5), g = list(mean = 0, sd = 1))
data1 <- rsimmix(n = 1000, unknownComp_weight = 0.6, list.comp, list.param)[['mixt.data']]
## Transform the known component of the admixture model into a Uniform(O,1) distribution:
list.comp <- list(f = NULL, g = 'norm')
list.param <- list(f = NULL, g = list(mean = 0, sd = 1))
data1_transfo <- knownComp_to_uniform(data = data1, comp.dist=list.comp, comp.param=list.param)
PatraSen_est_mix_model(data = data1_transfo, method = 'fixed',
c.n = 0.1*log(log(length(data1_transfo))), gridsize = 1000)$alp.hat
```
## The two-sample case
In the two-sample setting, the idea is to use the Inversion - Best Matching (IBM) approach. The IBM method sensures asymptotically normal estimators of the unknown quantities, which will be very useful in a testing perspective. However, it is important to note that such estimators are mostly biased when $F_1 \neq F_2$, and general one-sample estimation strategies such as [@Patra] or [@BV10] may be preferred to estimate the unknown component proportion in general settings (despite that this is much more time-consuming). In the latter case, one performs twice the estimation method, on each of the two samples under study.
### Under the null hypothesis $H_0: F_1 = F_2$
When we are under the null, @MilhaudIBM21 show that the estimatorsis consistent towards the true parameter values.
```{r}
## Simulate data:
list.comp <- list(f1 = 'norm', g1 = 'norm',
f2 = 'norm', g2 = 'norm')
list.param <- list(f1 = list(mean = 3, sd = 0.5), g1 = list(mean = 0, sd = 1),
f2 = list(mean = 3, sd = 0.5), g2 = list(mean = 5, sd = 2))
sample1 <- rsimmix(n=1700, unknownComp_weight=0.5, comp.dist = list(list.comp$f1,list.comp$g1),
comp.param=list(list.param$f1,list.param$g1))
sample2 <- rsimmix(n=1500, unknownComp_weight=0.7, comp.dist = list(list.comp$f2,list.comp$g2),
comp.param=list(list.param$f2,list.param$g2))
##### On a real-life example (unknown component densities, unknown mixture weights).
list.comp <- list(f1 = NULL, g1 = 'norm',
f2 = NULL, g2 = 'norm')
list.param <- list(f1 = NULL, g1 = list(mean = 0, sd = 1),
f2 = NULL, g2 = list(mean = 5, sd = 2))
## Estimate the mixture weights of the two admixture models (provide only hat(theta)_n):
estim <- IBM_estimProp(sample1 = sample1[['mixt.data']], sample2 = sample2[['mixt.data']],
known.prop = NULL, comp.dist = list.comp, comp.param = list.param,
with.correction = FALSE, n.integ = 1000)
estim[['prop.estim']]
```
Indeed, the two unknown proportions were consistently estimated.
### Under the alternative hypothesis $H_1: F_1 \neq F_2$
Estimators are consistent under $H_1$, although they can be (strongly) biased as compared to their true values.
```{r}
## Simulate data:
list.comp <- list(f1 = 'norm', g1 = 'norm',
f2 = 'norm', g2 = 'norm')
list.param <- list(f1 = list(mean = 1, sd = 0.5), g1 = list(mean = 0, sd = 1),
f2 = list(mean = 3, sd = 0.5), g2 = list(mean = 5, sd = 2))
sample1 <- rsimmix(n=1700, unknownComp_weight=0.5, comp.dist = list(list.comp$f1,list.comp$g1),
comp.param=list(list.param$f1,list.param$g1))
sample2 <- rsimmix(n=1500, unknownComp_weight=0.7, comp.dist = list(list.comp$f2,list.comp$g2),
comp.param=list(list.param$f2,list.param$g2))
##### On a real-life example (unknown component densities, unknown mixture weights).
list.comp <- list(f1 = NULL, g1 = 'norm',
f2 = NULL, g2 = 'norm')
list.param <- list(f1 = NULL, g1 = list(mean = 0, sd = 1),
f2 = NULL, g2 = list(mean = 5, sd = 2))
## Estimate the mixture weights of the two admixture models (provide only hat(theta)_n):
estim <- IBM_estimProp(sample1 = sample1[['mixt.data']], sample2 = sample2[['mixt.data']],
known.prop = NULL, comp.dist = list.comp, comp.param = list.param,
with.correction = FALSE, n.integ = 1000)
estim[['prop.estim']]
```
# Estimation of the unknown cumulative distribution function
Concerning the unknown cdf $F$, one usually estimate it thanks to the inversion formula
$$F(x) = \dfrac{L(x) - (1-p)G(x)}{p},$$
once $p$ has been consistenly estimated.
This is what is commonly called the decontaminated density of the unknown component. In the following, we propose to compare the two decontaminated densities obtained once the unknown quantities have been consistently estimated by the IBM approach. Note that we are under the null ($F_1=F_2$), and thus that the decontaminated densities should look similar.
```{r}
## Simulate data:
list.comp <- list(f1 = 'norm', g1 = 'norm',
f2 = 'norm', g2 = 'norm')
list.param <- list(f1 = list(mean = 3, sd = 0.5), g1 = list(mean = 0, sd = 1),
f2 = list(mean = 3, sd = 0.5), g2 = list(mean = 5, sd = 2))
sample1 <- rsimmix(n=1700, unknownComp_weight=0.5, comp.dist = list(list.comp$f1,list.comp$g1),
comp.param=list(list.param$f1,list.param$g1))
sample2 <- rsimmix(n=1500, unknownComp_weight=0.7, comp.dist = list(list.comp$f2,list.comp$g2),
comp.param=list(list.param$f2,list.param$g2))
## Estimate the mixture weight in each of the sample in real-life setting:
list.comp <- list(f1 = NULL, g1 = 'norm',
f2 = NULL, g2 = 'norm')
list.param <- list(f1 = NULL, g1 = list(mean = 0, sd = 1),
f2 = NULL, g2 = list(mean = 5, sd = 2))
estimate <- IBM_estimProp(sample1[['mixt.data']], sample2[['mixt.data']], comp.dist = list.comp,
comp.param = list.param, with.correction = FALSE, n.integ = 1000)
## Determine the decontaminated version of the unknown density by inversion:
res1 <- decontaminated_density(sample1 = sample1[['mixt.data']], comp.dist = list.comp[1:2],
comp.param = list.param[1:2], estim.p = estimate$prop.estim[1])
res2 <- decontaminated_density(sample1 = sample2[['mixt.data']], comp.dist = list.comp[3:4],
comp.param = list.param[3:4], estim.p = estimate$prop.estim[2])
plot(x = res1, type = "l", x_val = seq(from=-1,to=6,length.out=40), add_plot = FALSE)
plot(x = res2, type="l", col="red", x_val = seq(from=-1,to=6,length.out=40), add_plot = TRUE)
```
# References