Bayesian Linear Mixed Models (Conventional)



In [20]:

    
using JWAS,JWAS.Datasets,DataFrames,CSV



In [2]:

    
phenofile  = Datasets.dataset("example","phenotypes.txt")
phenotypes = CSV.read(phenofile,delim = ',',header=true,missingstrings=["NA"]);



In [3]:

    
first(phenotypes,5)









    Out[3]:




5 rows × 9 columns
ID y1 y2 y3 x1 x2 x3 dam weights
String Float64⍰ Float64⍰ Float64⍰ Float64 Float64 String String⍰ Float64
1 a1 -0.06 3.58 -1.18 0.9 2.0 m missing 1.0
2 a3 -2.07 3.19 1.0 0.7 2.0 f missing 1.0
3 a4 -2.63 6.97 -0.83 0.6 1.0 m a2 1.0
4 a5 2.31 missing -1.52 0.4 2.0 m a2 1.0
5 a7 missing missing missing 0.1 0.1 m a3 1.0

Univariate Linear Mixed Models (Conventional)



In [4]:

    
model_equation1  ="y1 = intercept + x1*x3 + x2 + x3";



In [5]:

    
model1 = build_model(model_equation1);



In [6]:

    
set_covariate(model1,"x1");



In [7]:

    
set_random(model1,"x2");



In [8]:

    
out1=runMCMC(model1,phenotypes);









    



Checking phenotypes...
Individual IDs (strings) are provided in the first column of the phenotypic data.
Phenotypes for all traits included in the model for individual a7 in the row 5 are missing. This record is deleted.
Prior information for residual variance is not provided and is generated from the data.
Prior information for random effect variance is not provided and is generated from the data.

A Linear Mixed Model was build using model equations:

y1 = intercept + x1*x3 + x2 + x3

Model Information:

Term            C/F          F/R            nLevels
intercept       factor       fixed                1
x1*x3           interaction  fixed                2
x2              factor       random               2
x3              factor       fixed                2

MCMC Information:

methods                   conventional (no markers)
chain_length                                    100
burnin                                            0
estimateScale                                 false
starting_value                                false
printout_frequency                              101
output_samples_frequency                          1
constraint                                    false
missing_phenotypes                             true
update_priors_frequency                           0
seed                                          false

Hyper-parameters Information:

random effect variances (y1:x2):              [1.253]
residual variances:                           1.253

Degree of freedom for hyper-parameters:

residual variances:                           4.000
random effect variances:                      5.000



The file MCMC_samples_residual_variance.txt already exists!!! It is overwritten by the new output.
The file MCMC_samples_y1.x2_variances.txt already exists!!! It is overwritten by the new output.






    



running MCMC for conventional (no markers)...100%|██████| Time: 0:00:01






    




The version of Julia and Platform in use:

Julia Version 1.3.1
Commit 2d5741174c (2019-12-30 21:36 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin18.6.0)
  CPU: Intel(R) Core(TM) i7-7700K CPU @ 4.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)


The analysis has finished. Results are saved in the returned variable and text files. MCMC samples are saved in text files.



In [9]:

    
keys(out1)









    Out[9]:





Base.KeySet for a Dict{Any,Any} with 2 entries. Keys:
  "location parameters"
  "residual variance"



In [10]:

    
out1["location parameters"]









    Out[10]:




7 rows × 5 columns
Trait Effect Level Estimate Std_Error
Any Any Any Any Any
1 y1 intercept intercept -0.68565 2.15387
2 y1 x1*x3 x1 * m -3.53323 2.7287
3 y1 x1*x3 x1 * f 2.51169 7.1064
4 y1 x2 2.0 1.01359 1.53674
5 y1 x2 1.0 -1.27899 1.43765
6 y1 x3 m 2.54268 3.66508
7 y1 x3 f -4.14117 3.3524



In [11]:

    
out1["residual variance"]









    Out[11]:




1 rows × 3 columns
Covariance Estimate Std_Error
Any Any Any
1 y1_y1 1.53495 1.65136

Multivariate Linear Mixed Models (Conventional)



In [12]:

    
model_equation2 ="y1 = intercept + x1 + x3
                  y2 = intercept + x1 + x2 + x3
                  y3 = intercept + x1 + x1*x3 + x2";



In [13]:

    
model2 = build_model(model_equation2);



In [14]:

    
set_covariate(model2,"x1");



In [15]:

    
set_random(model2,"x2");









    



x2 is not found in model equation 1.



In [16]:

    
out2=runMCMC(model2,phenotypes,chain_length=5000,output_samples_frequency=100);









    



Checking phenotypes...
Individual IDs (strings) are provided in the first column of the phenotypic data.
Phenotypes for all traits included in the model for individual a7 in the row 5 are missing. This record is deleted.
Prior information for residual variance is not provided and is generated from the data.
Prior information for random effect variance is not provided and is generated from the data.

A Linear Mixed Model was build using model equations:

y1 = intercept + x1 + x3
y2 = intercept + x1 + x2 + x3
y3 = intercept + x1 + x1*x3 + x2

Model Information:

Term            C/F          F/R            nLevels
intercept       factor       fixed                1
x1              covariate    fixed                1
x3              factor       fixed                2
x2              factor       random               2
x1*x3           interaction  fixed                2

MCMC Information:

methods                   conventional (no markers)
chain_length                                   5000
burnin                                            0
estimateScale                                 false
starting_value                                false
printout_frequency                             5001
output_samples_frequency                        100
constraint                                    false
missing_phenotypes                             true
update_priors_frequency                           0
seed                                          false

Hyper-parameters Information:

random effect variances (y2:x2,y3:x2):
 1.081  0.0  
 0.0    0.316
residual variances:           
 1.253  0.0    0.0  
 0.0    1.081  0.0  
 0.0    0.0    0.316

Degree of freedom for hyper-parameters:

residual variances:                           7.000
random effect variances:                      6.000



The file MCMC_samples_residual_variance.txt already exists!!! It is overwritten by the new output.
The file MCMC_samples_y2.x2_y3.x2_variances.txt already exists!!! It is overwritten by the new output.






    



running MCMC for conventional (no markers)...100%|██████| Time: 0:00:03






    




The version of Julia and Platform in use:

Julia Version 1.3.1
Commit 2d5741174c (2019-12-30 21:36 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin18.6.0)
  CPU: Intel(R) Core(TM) i7-7700K CPU @ 4.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)


The analysis has finished. Results are saved in the returned variable and text files. MCMC samples are saved in text files.



In [17]:

    
keys(out2)









    Out[17]:





Base.KeySet for a Dict{Any,Any} with 2 entries. Keys:
  "location parameters"
  "residual variance"



In [18]:

    
out2["location parameters"]









    Out[18]:




16 rows × 5 columns
Trait Effect Level Estimate Std_Error
Any Any Any Any Any
1 y1 intercept intercept 58.4043 34.3336
2 y1 x1 x1 -4.8379 4.53666
3 y1 x3 m -55.5406 34.0001
4 y1 x3 f -57.0454 34.0681
5 y2 intercept intercept 15.2137 14.6064
6 y2 x1 x1 -11.9029 8.04085
7 y2 x2 2.0 0.218078 1.02655
8 y2 x2 1.0 0.142865 1.07977
9 y2 x3 m -1.28758 15.4774
10 y2 x3 f -4.10139 15.3815
11 y3 intercept intercept -1.66267 1.01476
12 y3 x1 x1 -15.347 14.4518
13 y3 x1*x3 x1 * m 16.1353 14.3822
14 y3 x1*x3 x1 * f 19.3425 14.3084
15 y3 x2 2.0 -0.119527 0.59711
16 y3 x2 1.0 0.143961 0.493966



In [19]:

    
out2["residual variance"]









    Out[19]:




9 rows × 3 columns
Covariance Estimate Std_Error
Any Any Any
1 y1_y1 3.41408 2.91564
2 y1_y2 -0.236925 2.33926
3 y1_y3 -0.175038 0.861852
4 y2_y1 -0.236925 2.33926
5 y2_y2 1.53161 1.56651
6 y2_y3 0.0333331 0.479879
7 y3_y1 -0.175038 0.861852
8 y3_y2 0.0333331 0.479879
9 y3_y3 0.303169 0.273024

	ID	y1	y2	y3	x1	x2	x3	dam	weights
	String	Float64⍰	Float64⍰	Float64⍰	Float64	Float64	String	String⍰	Float64
1	a1	-0.06	3.58	-1.18	0.9	2.0	m	missing	1.0
2	a3	-2.07	3.19	1.0	0.7	2.0	f	missing	1.0
3	a4	-2.63	6.97	-0.83	0.6	1.0	m	a2	1.0
4	a5	2.31	missing	-1.52	0.4	2.0	m	a2	1.0
5	a7	missing	missing	missing	0.1	0.1	m	a3	1.0

	Trait	Effect	Level	Estimate	Std_Error
	Any	Any	Any	Any	Any
1	y1	intercept	intercept	-0.68565	2.15387
2	y1	x1*x3	x1 * m	-3.53323	2.7287
3	y1	x1*x3	x1 * f	2.51169	7.1064
4	y1	x2	2.0	1.01359	1.53674
5	y1	x2	1.0	-1.27899	1.43765
6	y1	x3	m	2.54268	3.66508
7	y1	x3	f	-4.14117	3.3524

	Trait	Effect	Level	Estimate	Std_Error
	Any	Any	Any	Any	Any
1	y1	intercept	intercept	58.4043	34.3336
2	y1	x1	x1	-4.8379	4.53666
3	y1	x3	m	-55.5406	34.0001
4	y1	x3	f	-57.0454	34.0681
5	y2	intercept	intercept	15.2137	14.6064
6	y2	x1	x1	-11.9029	8.04085
7	y2	x2	2.0	0.218078	1.02655
8	y2	x2	1.0	0.142865	1.07977
9	y2	x3	m	-1.28758	15.4774
10	y2	x3	f	-4.10139	15.3815
11	y3	intercept	intercept	-1.66267	1.01476
12	y3	x1	x1	-15.347	14.4518
13	y3	x1*x3	x1 * m	16.1353	14.3822
14	y3	x1*x3	x1 * f	19.3425	14.3084
15	y3	x2	2.0	-0.119527	0.59711
16	y3	x2	1.0	0.143961	0.493966

	Covariance	Estimate	Std_Error
	Any	Any	Any
1	y1_y1	3.41408	2.91564
2	y1_y2	-0.236925	2.33926
3	y1_y3	-0.175038	0.861852
4	y2_y1	-0.236925	2.33926
5	y2_y2	1.53161	1.56651
6	y2_y3	0.0333331	0.479879
7	y3_y1	-0.175038	0.861852
8	y3_y2	0.0333331	0.479879
9	y3_y3	0.303169	0.273024