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_equation  = "nwn = intercept +parity + parity*site + yr + geneticCode + age"


    

    
Out[3]:





"nwn = intercept +parity + parity*site + yr + geneticCode + age"

In [4]:

    

residual_variance = 2.97
model             = build_model(model_equation,residual_variance)

set_covariate(model,"age")

sow_variance      = 0.26
set_random(model,"parity",sow_variance);


    

In [5]:

    

outputMCMCsamples(model,"parity","age");


    

In [6]:

    

out=runMCMC(model,data,chain_length=50000,output_samples_frequency=100);


    

    




MCMC Information:
methods                        conventional (no markers)
chain_length                                  50000
estimatePi                                    false
constraint                                    false
missing_phenotypes                            false
starting_value                                false
output_samples_frequency                        100
printout_frequency                            50001
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:01

In [7]:

    

keys(out)


    

    
Out[7]:





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

In [8]:

    

out["Posterior mean of location parameters"]


    

    
Out[8]:





26x2 Array{Any,2}:
 "1:intercept : intercept"       85.074      
 "1:parity : 1"                   0.00084536 
 "1:parity : 2"                  -0.00194611 
 "1:parity : 3"                  -0.00156681 
 "1:parity : 4"                   0.000890018
 "1:parity : 5"                  -0.000997834
 "1:parity*site : 1 * 113"      294.825      
 "1:parity*site : 2 * 113"      131.797      
 "1:parity*site : 1 * 5"       -267.23       
 "1:parity*site : 1 * 13"      -341.117      
 "1:parity*site : 2 * 13"      -346.174      
 "1:parity*site : 3 * 13"      -351.203      
 "1:parity*site : 4 * 13"      -304.59       
 "1:parity*site : 5 * 13"      -329.748      
 "1:parity*site : 1 * 6"       -151.488      
 "1:yr : 2005"                   95.3079     
 "1:yr : 2006"                  262.351      
 "1:yr : 2008"                  -20.4951     
 "1:yr : 2009"                  -26.5392     
 "1:yr : 2004"                  135.895      
 "1:yr : 2012"                   41.6875     
 "1:geneticCode : PIC 1"       -521.642      
 "1:geneticCode : PIC 2"        167.298      
 "1:geneticCode : Commercial"    80.7843     
 "1:geneticCode : PIC C27"      -14.667      
 "1:age : age"                    3.02332    

In [9]:

    

oursample=out["MCMC samples for: 1:age"]


    

    
Out[9]:





1x501 Array{Any,2}:
 "1:age : age"  -0.348487  -0.060165  …  1.8925  2.3202  2.14014  2.02161

In [11]:

    

using Plots
sample4age = vec(oursample[2:end])
histogram(sample4age)


    

    
Out[11]:

In [12]:

    

out["MCMC samples for: variance of 1:parity"]


    

    
Out[12]:





1x500 Array{Float64,2}:
 0.0899984  0.0888872  0.0351931  …  0.256085  0.0780606  0.0828524