To illustrate the method I will use some data on the\u00a0concentration of Immunoglobin G (IgG) measured in a sample of 298 children aged under 6 years. IgG is a component of the immune system and can be found in the blood and other body fluids. It binds to the surfaces of bacteria\u00a0and viruses, immobilizing them and enabling phagocytic cells to recognise them. The normal range in adults is between 7 and 15 g\/L but in children the normal range is lower and varies with age. The data were used by Nuofaily and Jones (Applied Statistics<\/em> 2013;62:723-740) to illustrate quantile regression based on a generalized gamma distribution.<\/p>\n The idea behind the mixture method is to approximate the two-dimensional distribution of points by a finite mixture of normal distributions, as in the following illustration.<\/p>\n We can draw a slice through the mixture at a chosen age and obtain a flexible approximation to the conditional distribution by combining the component distributions. Alternatively we can find the conditional expectation, E(y|x), for any chosen age, x, by combining the conditional means of each of the components.<\/p>\n The fitted trend line is just the average of these expectations weighted by the probability of a point falling in that component and the probability of the particular x for that component. Repeating this calculation for lots of x’s\u00a0enables us to draw the smooth trend line, shown in blue in the illustration.<\/p>\n In the code below the priors and initial values are set and then the Bayesian mixture fitting program fits the model. The priors are set using the local r to denote the number of components in the mixture.<\/p>\n Now we can calculate the expectations and average them across the iterations. The chain of 5000 updates is thinned by taking every 10th from iteration 5 onwards in order to save time in the calculation.<\/p>\n Because we have a simulated Bayesian model fit it is a simple matter to find the average trend line and the 10% and 90% quantiles shown in red and blue in the plot below.<\/p>\n The choice of the number of components and setting of the priors together control the degree of smoothing and care should be taken when chosing these values. Perhaps by trying out different choices on simulated data to gain a feel for the degree of smoothing produced by a given set of priors.<\/p>\n In the paper by Noufaily and Jones the authors set themselves the slightly different problem of estimating the quantiles of the sample rather than smoothing the scatter plot. The fitted mixture model can be used to address this problem as well.<\/p>\n Imagine a vertical slice though the fitted bivariate components corresponding to some value x, we would obtain something similar to the picture below. The sum of these normal distributions produces represents the fitted conditional distribution for a particular age. The quantile that excludes q% will have q% of the sum of the areas in the blue region. All we need to do is to try a range of trial values of IgG until we find the value that excludes the required amount.<\/p>\n![\"qfig1\"](\"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/files\/2014\/09\/qfig1-300x218.png\")
<\/a><\/p>\n![\"qfig2\"](\"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/files\/2014\/09\/qfig2-300x218.png\")
<\/a><\/p>\n local r = 20<\/span>\r\n matrix R = J(2,2*`r',0)<\/span>\r\n forvalues j=1\/`r' {<\/span>\r\n matrix R[1,2*`j'-1] = 4<\/span>\r\n matrix R[2,2*`j'] = 4<\/span>\r\n }<\/span>\r\n matrix df = J(1,`r',16)<\/span>\r\n matrix M = J(2,`r',4)<\/span>\r\n matrix B = J(2,2*`r',0)<\/span>\r\n forvalues j=1\/`r' {<\/span>\r\n matrix B[1,2*`j'-1] = 0.1<\/span>\r\n matrix B[2,2*`j'] = 0.1<\/span>\r\n }<\/span>\r\n matrix ALPHA = J(1,`r',40\/`r')<\/span>\r\n matrix T = J(2,2*`r',0)<\/span>\r\n forvalues j=1\/`r' {<\/span>\r\n matrix T[1,2*`j'-1] = 1<\/span>\r\n matrix T[2,2*`j'] = 1<\/span>\r\n }<\/span>\r\n matrix C = J(298,1,1)<\/span>\r\n forvalues i=1\/298 {<\/span>\r\n matrix C[`i',1] = 1 + int(20*runiform())<\/span>\r\n }<\/span>\r\n matrix P = J(1,`r',1\/`r')<\/span>\r\n matrix MU = M\r\n <\/span>mixmnormal igg age, mu(MU) t(T) p(P) n(`r') m(M) b(B) r(R) \/\/\/<\/span>\r\n df(df) c(C) alpha(ALPHA) \/\/\/<\/span>\r\n burnin(1000) updates(5000) filename(\"igg.csv\")<\/span><\/pre>\n local r = 20<\/span>\r\n insheet using igg.csv, clear<\/span>\r\n merge 1:1 _n using igg.dta<\/span>\r\n range x 0 6 1001<\/span>\r\n foreach iter of numlist 5(10)5000 {<\/span>\r\n qui gen e`iter' = 0 in 1\/1001<\/span>\r\n qui gen w = 0 in 1\/1001<\/span>\r\n forvalues s=1\/`r' {<\/span>\r\n matrix T = J(2,2,0)<\/span>\r\n matrix T[1,1] = t`s'_1_1[`iter']<\/span>\r\n matrix T[2,1] = t`s'_2_1[`iter']<\/span>\r\n matrix T[1,2] = t`s'_1_2[`iter']<\/span>\r\n matrix T[2,2] = t`s'_2_2[`iter']<\/span>\r\n matrix V = syminv(T)<\/span>\r\n local sy = sqrt(V[1,1])<\/span>\r\n local sx = sqrt(V[2,2])<\/span>\r\n local rho = V[1,2]\/(`sy'*`sx')<\/span>\r\n local my = mu`s'_1[`iter']<\/span>\r\n local mx = mu`s'_2[`iter']<\/span>\r\n local p = p`s'[`iter']<\/span>\r\n qui replace e`iter' = e`iter' + `p'*(`my'+`rho'*`sy'* \/\/\/<\/span>\r\n (x - `mx')\/`sx')*normalden((x - `mx')\/`sx')<\/span>\r\n qui replace w = w + `p'*normalden((x - `mx')\/`sx')<\/span>\r\n }<\/span>\r\n qui replace e`iter' = e`iter'\/w<\/span>\r\n drop w<\/span>\r\n }<\/span>\r\n egen fmean = rowmean(e*)<\/span>\r\n egen fp10 = rowpctile(e*) , p(10)<\/span>\r\n gen fp10 = rowpctile(e*) , p(90)<\/span>\r\n twoway (scatter igg age) (line fmean x, lpat(solid) lcol(red)) \/\/\/<\/span>\r\n (line fp10 x, lpat(dash) lcol(blue)) (line fp90 x, lpat(dash) \/\/\/<\/span>\r\n lcol(blue)), leg(off) xtitle(Age) ytitle(IgG)<\/span><\/pre>\n![\"qfig3\"](\"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/files\/2014\/09\/qfig3-300x218.png\")
<\/a><\/p>\n
![\"qfig4\"](\"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/files\/2014\/09\/qfig4-300x218.png\")
<\/a><\/p>\n![\"qfig5\"](\"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/files\/2014\/09\/qfig5-300x218.png\")
<\/a><\/p>\n