In [1]:

    

using DataFrames,JWAS,JWAS.Datasets


    

In [2]:

    

phenofile = Datasets.dataset("testMME","data.txt")
data      = readtable(phenofile,separator = ',',header=true)


    

    
Out[2]:




sow
site
yr
age
geneticCode
parity
nwn
SYS
bw
1
100-113
113
2005
18
PIC 1
1
8
113_2005_WNTR
9.0
2
100-113
113
2006
18
PIC 1
2
12
113_2006_SPNG
8.0
3
100-5
5
2008
15
PIC 2
1
10
5_2008_ATMN
7.5
4
1000-5
5
2009
17
PIC 2
1
10
5_2009_SPNG
8.3
5
10000-13
13
2004
16
Commercial
1
9
13_2004_WNTR
4.3
6
10000-13
13
2004
18
Commercial
2
10
13_2004_SMMR
2.8
7
10000-13
13
2004
20
Commercial
3
11
13_2004_ATMN
3.9
8
10000-13
13
2005
18
Commercial
4
11
13_2005_SPNG
10.0
9
10000-13
13
2005
25
Commercial
5
7
13_2005_ATMN
4.0
10
10000-6
6
2012
17
PIC C27
1
12
6_2012_ATMN
8.9

In [3]:

    

model_equations    = "nwn = intercept+parity*site+ age
                      bw = intercept+parity*site+ age";

residual_variance = [10.0 2.0
                      2.0 1.0];


    

In [4]:

    

model = build_model(model_equations,residual_variance);


    

In [5]:

    

set_covariate(model,"age");


    

In [6]:

    

outputMCMCsamples(model,"age");
out=runMCMC(model,data,chain_length=10000,output_samples_frequency=100);


    

    




MCMC Information:
methods                        conventional (no markers)
chain_length                                  10000
estimatePi                                    false
constraint                                    false
missing_phenotypes                            false
starting_value                                false
output_samples_frequency                        100
printout_frequency                            10001
update_priors_frequency                           0

Degree of freedom for hyper-parameters:
residual variances:                           4.000
iid random effect variances:                  4.000
polygenic effect variances:                   4.000
marker effect variances:                      4.000



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

In [7]:

    

keys(out)


    

    
Out[7]:





Base.KeyIterator for a Dict{Any,Any} with 5 entries. Keys:
  "Posterior mean of residual covariance matrix"
  "MCMC samples for: 1:age"
  "MCMC samples for: 2:age"
  "MCMC samples for residual covariance matrix"
  "Posterior mean of location parameters"

In [8]:

    

out["Posterior mean of location parameters"]


    

    
Out[8]:





22x2 Array{Any,2}:
 "1:intercept : intercept"   59.801    
 "1:parity*site : 1 * 113"  -46.1203   
 "1:parity*site : 2 * 113"  -41.9941   
 "1:parity*site : 1 * 5"    -44.6901   
 "1:parity*site : 1 * 13"   -45.6208   
 "1:parity*site : 2 * 13"   -44.0001   
 "1:parity*site : 3 * 13"   -42.3559   
 "1:parity*site : 4 * 13"   -42.8851   
 "1:parity*site : 5 * 13"   -44.7483   
 "1:parity*site : 1 * 6"    -42.3547   
 "1:age : age"               -0.322195 
 "2:intercept : intercept"  -10.4375   
 "2:parity*site : 1 * 113"   19.9762   
 "2:parity*site : 2 * 113"   19.0112   
 "2:parity*site : 1 * 5"     18.8183   
 "2:parity*site : 1 * 13"    15.241    
 "2:parity*site : 2 * 13"    13.7963   
 "2:parity*site : 3 * 13"    14.9657   
 "2:parity*site : 4 * 13"    21.0389   
 "2:parity*site : 5 * 13"    15.217    
 "2:parity*site : 1 * 6"     19.8494   
 "2:age : age"               -0.0310703

In [9]:

    

outsample=out["MCMC samples for: 2:age"]


    

    
Out[9]:





1x101 Array{Any,2}:
 "2:age : age"  0.0406323  -0.419902  …  -0.162842  -0.153645  -0.0193681

In [ ]:

    

using Plots
histogram(outsample[2:end])


    

In [ ]: