Load packages
library("MCMCglmm")
library("litterDiallel")
library("PLMcctools")
library("xtable")
Load data, add column
data("litters")
litters$inbred <- ifelse(litters$Dam_Founder == litters$Sire_Founder, 1, 0)
Generate design matrices
matrices <- diallelMatrixMaker(data=litters, dam.col.name="Dam_Founder", sire.col.name="Sire_Founder",
batch.col.name="YearMonth", batch.1.col.name="litterorder")
M.matrices <- diallelMatrixMakeAndRotate(litters, dam.col.name="Dam_Founder",
sire.col.name="Sire_Founder", batch.col.name="YearMonth",
batch.1.col.name="litterorder")[[1]]
t.matrices <- diallelMatrixMakeAndRotate(litters, dam.col.name="Dam_Founder",
sire.col.name="Sire_Founder", batch.col.name="YearMonth",
batch.1.col.name="litterorder")[[2]]
Define design matrices for each random effects class
add.mat <- matrices$add.mat
mat.mat <- matrices$mat.mat
dam.mat <- matrices$dam.mat
sire.mat <- matrices$sire.mat
inbred.mat <- matrices$inbred.mat
jk.mat <- matrices$jk.mat
asymm.mat <- matrices$asymm.mat
batch.mat <- matrices$batch.mat
batch.1.mat <- matrices$batch.1.mat
Define reduced-dimension design matrices to pre-center variable estimates in each class
See: Lenarcic et al., 2012, Genetics; and Crowley et al., 2014, Genetics
t.add.mat <- t.matrices$t.add.mat
t.mat.mat <- t.matrices$t.mat.mat
t.dam.mat <- t.matrices$t.dam.mat
t.sire.mat <- t.matrices$t.sire.mat
t.inbred.mat <- t.matrices$t.inbred.mat
t.jk.mat <- t.matrices$t.jk.mat
t.asymm.mat <- t.matrices$t.asymm.mat
t.batch.mat <- t.matrices$t.batch.mat
t.batch.1.mat <- t.matrices$t.batch.1.mat
Define priors
priors <- list(B=list(V=1e+03*diag(3), mu=c(0,0,0)),
R = list(V = diag(1), nu = 0.002),
G = list(
G1=list(V=1, nu=0.002),
G2=list(V=1, nu=0.002),
G3=list(V=1, nu=0.002),
G4=list(V=1, nu=0.002),
G5=list(V=1, nu=0.002),
G6=list(V=1, nu=0.002),
G7=list(V=1, nu=0.002)))
Fit the model
fit1 <- MCMCglmm(SqrtWeaned ~ 1 + inbred + litternum,
random = ~ idv(t.add.mat) + idv(t.mat.mat) + idv(t.inbred.mat) +
idv(t.jk.mat) + idv(t.asymm.mat) +
idv(t.batch.mat) + idv(t.batch.1.mat),
prior=priors,
nitt=nitt, burnin = burnin,
thin=thin, pr=TRUE,
family = "gaussian", data=litters)
#> Warning: 'cBind' is deprecated.
#> Since R version 3.2.0, base's cbind() should work fine with S4 objects
#>
#> MCMC iteration = 0
#>
#> MCMC iteration = 1000
Retrieve MCMC chains
allChains <- as.mcmc(cbind(fit1$Sol,fit1$VCV))
allChains.fixed <- allChains[,c("(Intercept)", "inbred", "litternum")]
colnames(allChains.fixed) <- c("mu", "inbred", "litternum")
allChains.add <- allChains[, grep(colnames(allChains), pattern="^t.add.mat[^.]", value=TRUE)]
allChains.mat <- allChains[, grep(colnames(allChains), pattern="^t.mat.mat[^.]", value=TRUE)]
allChains.inbred <- allChains[, grep(colnames(allChains), pattern="^t.inbred.mat[^.]", value=TRUE)]
allChains.jk <- allChains[, grep(colnames(allChains), pattern="^t.jk.mat[^.]", value=TRUE)]
allChains.asymm <- allChains[, grep(colnames(allChains), pattern="^t.asymm.mat[^.]", value=TRUE)]
allChains.batch <- allChains[, grep(colnames(allChains), pattern="^t.batch.mat[^.]", value=TRUE)]
allChains.batch.1 <- allChains[, grep(colnames(allChains), pattern="^t.batch.1.mat[^.]", value=TRUE)]
allChains.sigma <- allChains[, c('t.add.mat.', 't.mat.mat.', 't.inbred.mat.',
't.jk.mat.', 't.asymm.mat.',
't.batch.mat.', 't.batch.1.mat.', 'units')]
Rotate parameter estimates back to full parameter space
M.add <- M.matrices$M.add
M.mat <- M.matrices$M.mat
M.inbred <- M.matrices$M.inbred
M.jk <- M.matrices$M.jk
M.asymm <- M.matrices$M.asymm
M.batch <- M.matrices$M.batch
M.batch.1 <- M.matrices$M.batch.1
allChains.add <- allChains.add %*% t(M.add)
allChains.mat <- allChains.mat %*% t(M.mat)
allChains.inbred <- allChains.inbred %*% t(M.inbred)
allChains.jk <- allChains.jk %*% t(M.jk)
allChains.asymm <- allChains.asymm %*% t(M.asymm)
allChains.batch <- allChains.batch %*% t(M.batch)
allChains.batch.1 <- allChains.batch.1 %*% t(M.batch.1)
colnames(allChains.add) <- colnames(add.mat)
colnames(allChains.mat) <- colnames(mat.mat)
colnames(allChains.inbred) <- colnames(inbred.mat)
colnames(allChains.jk) <- colnames(jk.mat)
colnames(allChains.asymm) <- colnames(asymm.mat)
colnames(allChains.batch) <- colnames(batch.mat)
colnames(allChains.batch.1) <- colnames(batch.1.mat)
Edit parameter class names, combine chains
colnames(allChains.sigma) <- c("tau-sq-add",
"tau-sq-maternal",
"tau-sq-inbred",
"tau-sq-jk",
"tau-sq-asymm",
"tau-sq-year-month",
"tau-sq-litterorder",
"sigma-sq")
allChains <- as.mcmc(cbind(allChains.fixed, allChains.add,
allChains.mat, allChains.inbred,
allChains.jk, allChains.asymm,
allChains.batch, allChains.batch.1,
allChains.sigma))
Plot MCMC traces
#plot(allChains[,1])
Plot diallel effect estimates
jk.ordering <- c('v:B6;AJ', 'v:129;AJ', 'v:NOD;AJ', 'v:NZO;AJ', 'v:CAST;AJ', 'v:PWK;AJ', 'v:WSB;AJ',
'v:129;B6', 'v:NOD;B6', 'v:NZO;B6', 'v:CAST;B6', 'v:PWK;B6', 'v:WSB;B6',
'v:NOD;129', 'v:NZO;129', 'v:CAST;129', 'v:PWK;129', 'v:WSB;129',
'v:NZO;NOD', 'v:CAST;NOD', 'v:PWK;NOD', 'v:WSB;NOD',
'v:CAST;NZO', 'v:PWK;NZO', 'v:WSB;NZO', 'v:PWK;CAST', 'v:WSB;CAST', 'v:WSB;PWK')
asymm.ordering <- c('w:B6;AJ', 'w:129;AJ', 'w:NOD;AJ', 'w:NZO;AJ', 'w:CAST;AJ', 'w:PWK;AJ', 'w:WSB;AJ',
'w:129;B6', 'w:NOD;B6', 'w:NZO;B6', 'w:CAST;B6', 'w:PWK;B6', 'w:WSB;B6',
'w:NOD;129', 'w:NZO;129', 'w:CAST;129', 'w:PWK;129', 'w:WSB;129',
'w:NZO;NOD', 'w:CAST;NOD', 'w:PWK;NOD', 'w:WSB;NOD',
'w:CAST;NZO', 'w:PWK;NZO', 'w:WSB;NZO', 'w:PWK;CAST', 'w:WSB;CAST', 'w:WSB;PWK')
varcolors <- c("#A6CEE3", "#B2DF8A", "#FB9A99",
"#FDBF6F", "#CAB2D6", "#D2B48C")
varcolors.y0 <- rev(c(0.25, 8.25, 9.25, 17.25))
varcolors.y1 <- rev(c(7.75, 8.75, 16.75, 24.75))
xlim=c(-0.3, 0.3);
var.labels.all <- varnames(allChains)
var.labels <- c(colnames(add.mat), colnames(mat.mat), "inbred", colnames(inbred.mat))
plot.hpd(allChains[,var.labels], main="general effects", xlim=xlim); abline(v=0, col="grey")
segments(x0=xlim[1]*0.85, x1=xlim[1]*0.85, y0=varcolors.y0, y1=varcolors.y1, col=varcolors[1:4], lwd=7, lend=1)
var.labels <- jk.ordering
plot.hpd(allChains[,var.labels], main="strainpair-specific", xlim=xlim); abline(v=0, col="grey")
segments(x0=xlim[1]*0.85, x1=xlim[1]*0.85, y0=0.25, y1=27.75, col=varcolors[5], lwd=7, lend=1)
var.labels <- asymm.ordering
plot.hpd(allChains[,var.labels], main="strainpair-specific", xlim=xlim); abline(v=0, col="grey")
segments(x0=xlim[1]*0.85, x1=xlim[1]*0.85, y0=0.25, y1=27.75, col=varcolors[6], lwd=7, lend=1)
Plot covariate effect estimates
var.labels.all <- varnames(allChains)
var.labels <- c(colnames(batch.mat))[1:24]
plot.hpd(allChains[,var.labels], main="wean date", xlim=xlim); abline(v=0, col="grey")
var.labels <- c(colnames(batch.mat))[25:48]
plot.hpd(allChains[,var.labels], main="wean date", xlim=xlim); abline(v=0, col="grey")
var.labels <- c("litternum", colnames(batch.1.mat))
plot.hpd(allChains[,var.labels], main="litter order", xlim=xlim); abline(v=0, col="grey")
Print tables
#hpd.table1 <- summary(allChains)[[1]]
hpd.table2 <- summary(allChains)[[2]]
#print(xtable(hpd.table1), type="html")
print(xtable(hpd.table2, digits=4), type="html")
|
|
2.5%
|
25%
|
50%
|
75%
|
97.5%
|
|
mu
|
2.4177
|
2.4676
|
2.4943
|
2.5205
|
2.5832
|
|
inbred
|
-0.2153
|
-0.1953
|
-0.1839
|
-0.1729
|
-0.1517
|
|
litternum
|
-0.0606
|
-0.0454
|
-0.0378
|
-0.0304
|
-0.0192
|
|
additive:AJ
|
-0.0874
|
-0.0548
|
-0.0386
|
-0.0216
|
0.0090
|
|
additive:B6
|
0.0778
|
0.1167
|
0.1366
|
0.1558
|
0.1902
|
|
additive:129
|
-0.0624
|
-0.0303
|
-0.0151
|
-0.0000
|
0.0278
|
|
additive:NOD
|
0.1293
|
0.1622
|
0.1784
|
0.1941
|
0.2231
|
|
additive:NZO
|
0.0387
|
0.0735
|
0.0925
|
0.1119
|
0.1502
|
|
additive:CAST
|
-0.1801
|
-0.1455
|
-0.1258
|
-0.1052
|
-0.0632
|
|
additive:PWK
|
-0.1004
|
-0.0698
|
-0.0547
|
-0.0401
|
-0.0117
|
|
additive:WSB
|
-0.2094
|
-0.1743
|
-0.1606
|
-0.1448
|
-0.1110
|
|
maternal:AJ
|
-0.0422
|
-0.0035
|
0.0167
|
0.0354
|
0.0772
|
|
maternal:B6
|
0.0846
|
0.1247
|
0.1453
|
0.1649
|
0.2036
|
|
maternal:129
|
-0.0489
|
-0.0058
|
0.0147
|
0.0343
|
0.0784
|
|
maternal:NOD
|
0.0256
|
0.0673
|
0.0888
|
0.1125
|
0.1570
|
|
maternal:NZO
|
-0.0491
|
-0.0072
|
0.0182
|
0.0421
|
0.0895
|
|
maternal:CAST
|
-0.0915
|
-0.0525
|
-0.0294
|
-0.0088
|
0.0370
|
|
maternal:PWK
|
-0.1694
|
-0.1286
|
-0.1090
|
-0.0885
|
-0.0472
|
|
maternal:WSB
|
-0.1886
|
-0.1532
|
-0.1346
|
-0.1149
|
-0.0789
|
|
inbred:AJ
|
-0.0614
|
-0.0046
|
0.0275
|
0.0600
|
0.1323
|
|
inbred:B6
|
-0.0649
|
-0.0015
|
0.0338
|
0.0714
|
0.1431
|
|
inbred:129
|
-0.2476
|
-0.1700
|
-0.1277
|
-0.0860
|
-0.0075
|
|
inbred:NOD
|
0.0071
|
0.0890
|
0.1358
|
0.1774
|
0.2545
|
|
inbred:NZO
|
-0.0587
|
0.0003
|
0.0357
|
0.0709
|
0.1373
|
|
inbred:CAST
|
-0.1898
|
-0.1062
|
-0.0625
|
-0.0190
|
0.0461
|
|
inbred:PWK
|
-0.1026
|
-0.0374
|
-0.0078
|
0.0236
|
0.0881
|
|
inbred:WSB
|
-0.1375
|
-0.0626
|
-0.0294
|
0.0014
|
0.0615
|
|
v:129;AJ
|
-0.0948
|
-0.0325
|
-0.0102
|
0.0107
|
0.0543
|
|
v:129;B6
|
-0.0203
|
0.0159
|
0.0377
|
0.0658
|
0.1291
|
|
v:B6;AJ
|
-0.0626
|
-0.0167
|
0.0053
|
0.0250
|
0.0858
|
|
v:CAST;129
|
-0.1224
|
-0.0511
|
-0.0242
|
-0.0020
|
0.0390
|
|
v:CAST;AJ
|
-0.0823
|
-0.0229
|
0.0006
|
0.0219
|
0.0704
|
|
v:CAST;B6
|
-0.0556
|
-0.0166
|
0.0045
|
0.0264
|
0.0702
|
|
v:CAST;NOD
|
-0.1012
|
-0.0368
|
-0.0112
|
0.0100
|
0.0525
|
|
v:CAST;NZO
|
-0.0976
|
-0.0301
|
-0.0061
|
0.0162
|
0.0703
|
|
v:NOD;129
|
-0.1511
|
-0.0706
|
-0.0411
|
-0.0163
|
0.0213
|
|
v:NOD;AJ
|
-0.0645
|
-0.0177
|
0.0039
|
0.0270
|
0.0795
|
|
v:NOD;B6
|
-0.0550
|
-0.0126
|
0.0088
|
0.0337
|
0.0936
|
|
v:NZO;129
|
-0.0328
|
0.0076
|
0.0317
|
0.0614
|
0.1378
|
|
v:NZO;AJ
|
-0.0644
|
-0.0209
|
0.0028
|
0.0249
|
0.0807
|
|
v:NZO;B6
|
-0.0987
|
-0.0338
|
-0.0101
|
0.0129
|
0.0608
|
|
v:NZO;NOD
|
-0.0520
|
-0.0100
|
0.0135
|
0.0385
|
0.1058
|
|
v:PWK;129
|
-0.0996
|
-0.0436
|
-0.0194
|
0.0017
|
0.0437
|
|
v:PWK;AJ
|
-0.0793
|
-0.0223
|
-0.0008
|
0.0198
|
0.0673
|
|
v:PWK;B6
|
-0.0579
|
-0.0112
|
0.0080
|
0.0315
|
0.0928
|
|
v:PWK;CAST
|
-0.0665
|
-0.0152
|
0.0056
|
0.0270
|
0.0705
|
|
v:PWK;NOD
|
-0.0624
|
-0.0168
|
0.0038
|
0.0272
|
0.0726
|
|
v:PWK;NZO
|
-0.0852
|
-0.0238
|
0.0009
|
0.0260
|
0.0893
|
|
v:WSB;129
|
-0.0528
|
-0.0075
|
0.0133
|
0.0366
|
0.0967
|
|
v:WSB;AJ
|
-0.0970
|
-0.0392
|
-0.0151
|
0.0067
|
0.0461
|
|
v:WSB;B6
|
-0.0483
|
-0.0108
|
0.0103
|
0.0346
|
0.0905
|
|
v:WSB;CAST
|
-0.1360
|
-0.0624
|
-0.0344
|
-0.0099
|
0.0302
|
|
v:WSB;NOD
|
-0.0395
|
0.0012
|
0.0238
|
0.0478
|
0.1074
|
|
v:WSB;NZO
|
-0.0618
|
-0.0173
|
0.0047
|
0.0289
|
0.0865
|
|
v:WSB;PWK
|
-0.0695
|
-0.0264
|
-0.0052
|
0.0160
|
0.0621
|
|
w:129;AJ
|
-0.1009
|
-0.0289
|
0.0103
|
0.0458
|
0.1297
|
|
w:129;B6
|
-0.0416
|
0.0206
|
0.0541
|
0.0898
|
0.1611
|
|
w:B6;AJ
|
-0.0774
|
-0.0069
|
0.0306
|
0.0675
|
0.1475
|
|
w:CAST;129
|
-0.0228
|
0.0521
|
0.0900
|
0.1303
|
0.2181
|
|
w:CAST;AJ
|
-0.0681
|
0.0005
|
0.0368
|
0.0763
|
0.1547
|
|
w:CAST;B6
|
-0.1620
|
-0.0858
|
-0.0513
|
-0.0155
|
0.0502
|
|
w:CAST;NOD
|
-0.1431
|
-0.0657
|
-0.0300
|
0.0098
|
0.0788
|
|
w:CAST;NZO
|
-0.1528
|
-0.0595
|
-0.0153
|
0.0280
|
0.1114
|
|
w:NOD;129
|
-0.1945
|
-0.1023
|
-0.0626
|
-0.0231
|
0.0463
|
|
w:NOD;AJ
|
-0.0753
|
-0.0076
|
0.0337
|
0.0727
|
0.1634
|
|
w:NOD;B6
|
-0.0879
|
-0.0123
|
0.0257
|
0.0625
|
0.1353
|
|
w:NZO;129
|
-0.1408
|
-0.0525
|
-0.0160
|
0.0260
|
0.1060
|
|
w:NZO;AJ
|
-0.1232
|
-0.0457
|
-0.0076
|
0.0317
|
0.1086
|
|
w:NZO;B6
|
-0.1151
|
-0.0375
|
0.0030
|
0.0418
|
0.1298
|
|
w:NZO;NOD
|
-0.1042
|
-0.0191
|
0.0208
|
0.0624
|
0.1425
|
|
w:PWK;129
|
-0.0779
|
-0.0085
|
0.0251
|
0.0610
|
0.1410
|
|
w:PWK;AJ
|
-0.1306
|
-0.0631
|
-0.0261
|
0.0072
|
0.0713
|
|
w:PWK;B6
|
-0.1637
|
-0.0894
|
-0.0561
|
-0.0228
|
0.0373
|
|
w:PWK;CAST
|
-0.0246
|
0.0406
|
0.0763
|
0.1130
|
0.1837
|
|
w:PWK;NOD
|
-0.2642
|
-0.1780
|
-0.1366
|
-0.1000
|
-0.0302
|
|
w:PWK;NZO
|
-0.1226
|
-0.0374
|
0.0070
|
0.0530
|
0.1395
|
|
w:WSB;129
|
-0.0359
|
0.0301
|
0.0669
|
0.1034
|
0.1767
|
|
w:WSB;AJ
|
-0.1120
|
-0.0431
|
-0.0105
|
0.0227
|
0.0905
|
|
w:WSB;B6
|
-0.0810
|
-0.0130
|
0.0221
|
0.0535
|
0.1179
|
|
w:WSB;CAST
|
-0.1802
|
-0.0990
|
-0.0600
|
-0.0256
|
0.0503
|
|
w:WSB;NOD
|
-0.0194
|
0.0587
|
0.0969
|
0.1340
|
0.2102
|
|
w:WSB;NZO
|
-0.1503
|
-0.0665
|
-0.0264
|
0.0120
|
0.0873
|
|
w:WSB;PWK
|
-0.2069
|
-0.1291
|
-0.0930
|
-0.0547
|
0.0115
|
|
YearMonth_2008-03
|
-0.0427
|
0.0049
|
0.0317
|
0.0593
|
0.1199
|
|
YearMonth_2008-04
|
-0.0368
|
0.0109
|
0.0367
|
0.0670
|
0.1359
|
|
YearMonth_2008-05
|
-0.0273
|
0.0145
|
0.0398
|
0.0659
|
0.1311
|
|
YearMonth_2008-06
|
-0.0425
|
0.0073
|
0.0320
|
0.0585
|
0.1229
|
|
YearMonth_2008-07
|
-0.0388
|
0.0110
|
0.0362
|
0.0657
|
0.1277
|
|
YearMonth_2008-08
|
-0.0740
|
-0.0201
|
0.0060
|
0.0300
|
0.0832
|
|
YearMonth_2008-09
|
-0.0803
|
-0.0278
|
-0.0043
|
0.0207
|
0.0705
|
|
YearMonth_2008-10
|
-0.0649
|
-0.0161
|
0.0100
|
0.0353
|
0.0871
|
|
YearMonth_2008-11
|
-0.0752
|
-0.0219
|
0.0025
|
0.0289
|
0.0803
|
|
YearMonth_2008-12
|
-0.0700
|
-0.0167
|
0.0086
|
0.0366
|
0.0942
|
|
YearMonth_2009-01
|
-0.0752
|
-0.0220
|
0.0030
|
0.0295
|
0.0838
|
|
YearMonth_2009-02
|
-0.0883
|
-0.0255
|
-0.0015
|
0.0220
|
0.0733
|
|
YearMonth_2009-03
|
-0.1003
|
-0.0441
|
-0.0159
|
0.0092
|
0.0636
|
|
YearMonth_2009-04
|
-0.0769
|
-0.0205
|
0.0053
|
0.0327
|
0.0896
|
|
YearMonth_2009-05
|
-0.0885
|
-0.0350
|
-0.0075
|
0.0173
|
0.0707
|
|
YearMonth_2009-06
|
-0.0762
|
-0.0247
|
0.0019
|
0.0259
|
0.0797
|
|
YearMonth_2009-07
|
-0.0579
|
-0.0079
|
0.0170
|
0.0424
|
0.0962
|
|
YearMonth_2009-08
|
-0.0588
|
-0.0163
|
0.0083
|
0.0320
|
0.0772
|
|
YearMonth_2009-09
|
-0.0263
|
0.0101
|
0.0304
|
0.0499
|
0.0870
|
|
YearMonth_2009-10
|
-0.0437
|
-0.0022
|
0.0187
|
0.0399
|
0.0801
|
|
YearMonth_2009-11
|
-0.0611
|
-0.0243
|
-0.0046
|
0.0138
|
0.0538
|
|
YearMonth_2009-12
|
-0.0070
|
0.0252
|
0.0420
|
0.0628
|
0.1027
|
|
YearMonth_2010-01
|
-0.0464
|
-0.0095
|
0.0110
|
0.0312
|
0.0705
|
|
YearMonth_2010-02
|
-0.0910
|
-0.0433
|
-0.0214
|
-0.0006
|
0.0419
|
|
YearMonth_2010-03
|
-0.0413
|
-0.0034
|
0.0157
|
0.0345
|
0.0729
|
|
YearMonth_2010-04
|
-0.0924
|
-0.0532
|
-0.0319
|
-0.0137
|
0.0196
|
|
YearMonth_2010-05
|
-0.1514
|
-0.1029
|
-0.0812
|
-0.0602
|
-0.0253
|
|
YearMonth_2010-06
|
-0.1018
|
-0.0661
|
-0.0467
|
-0.0291
|
0.0057
|
|
YearMonth_2010-07
|
-0.0632
|
-0.0271
|
-0.0063
|
0.0112
|
0.0508
|
|
YearMonth_2010-08
|
-0.1162
|
-0.0639
|
-0.0414
|
-0.0203
|
0.0207
|
|
YearMonth_2010-09
|
-0.0273
|
0.0080
|
0.0285
|
0.0478
|
0.0922
|
|
YearMonth_2010-10
|
-0.0926
|
-0.0479
|
-0.0247
|
-0.0019
|
0.0395
|
|
YearMonth_2010-11
|
-0.0968
|
-0.0527
|
-0.0284
|
-0.0057
|
0.0324
|
|
YearMonth_2010-12
|
-0.0952
|
-0.0501
|
-0.0273
|
-0.0063
|
0.0313
|
|
YearMonth_2011-01
|
-0.0687
|
-0.0283
|
-0.0046
|
0.0183
|
0.0545
|
|
YearMonth_2011-02
|
-0.0787
|
-0.0342
|
-0.0120
|
0.0096
|
0.0507
|
|
YearMonth_2011-03
|
-0.0653
|
-0.0280
|
-0.0081
|
0.0098
|
0.0431
|
|
YearMonth_2011-04
|
-0.0761
|
-0.0360
|
-0.0155
|
0.0047
|
0.0418
|
|
YearMonth_2011-05
|
-0.0658
|
-0.0281
|
-0.0081
|
0.0113
|
0.0526
|
|
YearMonth_2011-06
|
-0.0559
|
-0.0187
|
0.0009
|
0.0206
|
0.0581
|
|
YearMonth_2011-07
|
-0.0635
|
-0.0217
|
-0.0003
|
0.0196
|
0.0631
|
|
YearMonth_2011-08
|
-0.0596
|
-0.0193
|
0.0002
|
0.0195
|
0.0614
|
|
YearMonth_2011-09
|
-0.0767
|
-0.0325
|
-0.0096
|
0.0102
|
0.0503
|
|
YearMonth_2011-10
|
-0.0856
|
-0.0331
|
-0.0090
|
0.0128
|
0.0563
|
|
YearMonth_2011-11
|
-0.0492
|
-0.0034
|
0.0185
|
0.0407
|
0.0949
|
|
YearMonth_2011-12
|
-0.0704
|
-0.0194
|
0.0062
|
0.0308
|
0.0815
|
|
YearMonth_2012-01
|
-0.0619
|
-0.0097
|
0.0136
|
0.0425
|
0.0960
|
|
YearMonth_2012-02
|
-0.0978
|
-0.0389
|
-0.0104
|
0.0154
|
0.0731
|
|
order01
|
-0.1929
|
-0.1251
|
-0.0953
|
-0.0716
|
-0.0245
|
|
order02
|
-0.0391
|
0.0084
|
0.0286
|
0.0500
|
0.0908
|
|
order03
|
-0.0215
|
0.0110
|
0.0284
|
0.0448
|
0.0803
|
|
order04
|
-0.0257
|
0.0071
|
0.0249
|
0.0413
|
0.0777
|
|
order05
|
-0.0194
|
0.0157
|
0.0341
|
0.0529
|
0.0925
|
|
order06
|
-0.0829
|
-0.0400
|
-0.0170
|
0.0047
|
0.0508
|
|
order07
|
-0.0612
|
-0.0089
|
0.0156
|
0.0439
|
0.1075
|
|
order08
|
-0.0410
|
0.0139
|
0.0459
|
0.0814
|
0.1893
|
|
order09
|
-0.1176
|
-0.0417
|
-0.0058
|
0.0273
|
0.0983
|
|
order10
|
-0.1501
|
-0.0621
|
-0.0238
|
0.0116
|
0.0811
|
|
order11
|
-0.1605
|
-0.0574
|
-0.0181
|
0.0204
|
0.0971
|
|
order12
|
-0.1403
|
-0.0482
|
-0.0108
|
0.0242
|
0.1067
|
|
tau-sq-add
|
0.0053
|
0.0100
|
0.0145
|
0.0220
|
0.0538
|
|
tau-sq-maternal
|
0.0032
|
0.0062
|
0.0091
|
0.0138
|
0.0369
|
|
tau-sq-inbred
|
0.0010
|
0.0046
|
0.0080
|
0.0136
|
0.0406
|
|
tau-sq-jk
|
0.0004
|
0.0008
|
0.0015
|
0.0026
|
0.0063
|
|
tau-sq-asymm
|
0.0022
|
0.0043
|
0.0060
|
0.0081
|
0.0151
|
|
tau-sq-year-month
|
0.0007
|
0.0013
|
0.0018
|
0.0024
|
0.0045
|
|
tau-sq-litterorder
|
0.0010
|
0.0020
|
0.0033
|
0.0053
|
0.0142
|
|
sigma-sq
|
0.2374
|
0.2440
|
0.2476
|
0.2512
|
0.2577
|
Calculate dam_strain and sire_strain effects
allChains.df <- as.data.frame(as.matrix(allChains))
allChains.df$`dam:AJ` <- allChains.df[,"additive:AJ"] + allChains.df[, "maternal:AJ"]
allChains.df$`dam:B6` <- allChains.df[,"additive:B6"] + allChains.df[, "maternal:B6"]
allChains.df$`dam:129` <- allChains.df[,"additive:129"] + allChains.df[, "maternal:129"]
allChains.df$`dam:NOD` <- allChains.df[,"additive:NOD"] + allChains.df[, "maternal:NOD"]
allChains.df$`dam:NZO` <- allChains.df[,"additive:NZO"] + allChains.df[, "maternal:NZO"]
allChains.df$`dam:CAST` <- allChains.df[,"additive:CAST"] + allChains.df[, "maternal:CAST"]
allChains.df$`dam:PWK` <- allChains.df[,"additive:PWK"] + allChains.df[, "maternal:PWK"]
allChains.df$`dam:WSB` <- allChains.df[,"additive:WSB"] + allChains.df[, "maternal:WSB"]
allChains.df$`sire:AJ` <- allChains.df[,"additive:AJ"] - allChains.df[, "maternal:AJ"]
allChains.df$`sire:B6` <- allChains.df[,"additive:B6"] - allChains.df[, "maternal:B6"]
allChains.df$`sire:129` <- allChains.df[,"additive:129"] - allChains.df[, "maternal:129"]
allChains.df$`sire:NOD` <- allChains.df[,"additive:NOD"] - allChains.df[, "maternal:NOD"]
allChains.df$`sire:NZO` <- allChains.df[,"additive:NZO"] - allChains.df[, "maternal:NZO"]
allChains.df$`sire:CAST` <- allChains.df[,"additive:CAST"] - allChains.df[, "maternal:CAST"]
allChains.df$`sire:PWK` <- allChains.df[,"additive:PWK"] - allChains.df[, "maternal:PWK"]
allChains.df$`sire:WSB` <- allChains.df[,"additive:WSB"] - allChains.df[, "maternal:WSB"]
allChains.dam <- allChains.df[, grep(names(allChains.df), pattern="^dam:", value=TRUE)]
allChains.sire <- allChains.df[, grep(names(allChains.df), pattern="^sire:", value=TRUE)]
allChains.bup <- allChains
allChains <- as.mcmc(cbind(allChains.fixed, allChains.add, allChains.mat,
allChains.dam, allChains.sire, allChains.inbred,
allChains.jk, allChains.asymm,
allChains.batch, allChains.batch.1,
allChains.sigma))
Plot dam_strain and _strain MCMC traces
#plot(allChains[,colnames(dam.mat)]); plot(allChains[,colnames(sire.mat)])
Plot dam_strain and sire_strain effects
xlim=c(-0.4, 0.4);
var.labels.all <- varnames(allChains)
var.labels.new <- c("dam:AJ", "dam:B6", "dam:129", "dam:NOD", "dam:NZO", "dam:CAST", "dam:PWK", "dam:WSB",
"sire:AJ", "sire:B6", "sire:129", "sire:NOD", "sire:NZO", "sire:CAST", "sire:PWK", "sire:WSB")
var.labels <- c(colnames(dam.mat), colnames(sire.mat))#, "inbred", colnames(inbred.mat))
plot.hpd(allChains[,var.labels[1:8]], xlim=xlim, main="dam strain", main.line=1.75); abline(v=0, col="grey")
plot.hpd(allChains[,var.labels[9:16]], xlim=xlim, main="sire strain", main.line=1.75); abline(v=0, col="grey")
Plot VarComps
allChains.sigma.names <- c("tau-sq-add",
"tau-sq-maternal",
"tau-sq-inbred",
"tau-sq-jk",
"tau-sq-asymm",
"sigma-sq")
allChains.sigma <- allChains[,allChains.sigma.names]
var.df <- cbind(as.data.frame(allChains[, allChains.sigma.names]), adjust)
names(var.df)[names(var.df)=="var1"] <- "adjust"
var.df$total <- var.df$`tau-sq-add` + var.df$`tau-sq-maternal` + var.df$`tau-sq-inbred` +
var.df$`tau-sq-jk` + var.df$`tau-sq-asymm` +
var.df$`sigma-sq` + var.df$`adjust`
varcomp <- NULL
varcomp$additive <- var.df$`tau-sq-add` / var.df$`total`
varcomp$parental.sex <- var.df$`tau-sq-maternal` / var.df$`total`
varcomp$inbred <- var.df$`tau-sq-inbred` / var.df$`total`
varcomp$symmetric.epistatic <- var.df$`tau-sq-jk` / var.df$`total`
varcomp$asymmetric.epistatic <- var.df$`tau-sq-asymm` / var.df$`total`
varcomp$total.explained <- (var.df$`total` - var.df$`sigma-sq` - var.df$`adjust`)/ var.df$`total`
varcomp$noise <- (var.df$`sigma-sq` + var.df$`adjust`)/var.df$`total`
varcomp <- as.data.frame(varcomp)
varcomp.table <- cbind(mean=colMeans(varcomp),
lower.bound=HPDinterval(as.mcmc(varcomp))[,1],
upper.bound=HPDinterval(as.mcmc(varcomp))[,2])
varcolors <- c("#A6CEE3", "#B2DF8A", #"#FB9A99",
"#FDBF6F", "#CAB2D6", "#D2B48C")
varcolors.y0 <- c(7.5, 6.5, 5.5, 4.5, 3.5)
varcolors.y1 <- c(6.5, 5.5, 4.5, 3.5, 2.5)
psq <- data.frame(varcomp.table)
par(mar = c(5.1, 13.5, 2.1, 12.1))
rows <- nrow(psq)
plot(y = c(0.5, rows + 0.5), x = c(-0.06, 1.01), yaxs = "i",
xaxs = "i", col = "white", ylab = "", xlab = "variance explained",
yaxt = "n")
segments(x0=-0.05, x1=-0.05, y0=varcolors.y0, y1=varcolors.y1, col=varcolors, lwd=15, lend=1)
abline(v = c(0, 0.2, 0.4, 0.6, 0.8, 1), lty = 2, col = "grey")
abline(h = 2.5, lwd = 1.5)
segments(x0 = psq$lower.bound, x1 = psq$upper.bound, y0 = c(rows:1),
y1 = c(rows:1), lwd = 4, lend=1)
axis(side = 2, at = c(rows:1), labels = rownames(psq),
las = 1)
rightlabels <- paste(as.character(sprintf("%.2f", round(100 *
psq$mean, digits = 2))), "% (", as.character(sprintf("%.2f",
round(100 * psq$lower.bound, digits = 2))), "%, ", as.character(sprintf("%.2f",
round(100 * psq$upper.bound, digits = 2))), "%)", sep = "")
axis(side = 4, at = c(rows:1), labels = rightlabels, tick = FALSE,
las = 1)
points(x = psq$mean, y = c(rows:1), pch = 16, col = "black",
cex = 1.5)
points(x = psq$mean, y = c(rows:1), pch = 16, col = "white",
cex = 0.5)
Table of results
print(xtable(psq, digits=4), type="html")
|
|
mean
|
lower.bound
|
upper.bound
|
|
additive
|
0.0596
|
0.0124
|
0.1349
|
|
parental.sex
|
0.0389
|
0.0077
|
0.0930
|
|
inbred
|
0.0363
|
0.0010
|
0.0953
|
|
symmetric.epistatic
|
0.0067
|
0.0008
|
0.0183
|
|
asymmetric.epistatic
|
0.0224
|
0.0059
|
0.0458
|
|
total.explained
|
0.1639
|
0.0832
|
0.2739
|
|
noise
|
0.8361
|
0.7261
|
0.9168
|
Plot VarPs
allChains.df <- data.frame(allChains, check.names=FALSE)
allChains.df$`dam:AJ` <- allChains.df[,"additive:AJ"] + allChains.df[, "maternal:AJ"]
allChains.df$`dam:B6` <- allChains.df[,"additive:B6"] + allChains.df[, "maternal:B6"]
allChains.df$`dam:129` <- allChains.df[,"additive:129"] + allChains.df[, "maternal:129"]
allChains.df$`dam:NOD` <- allChains.df[,"additive:NOD"] + allChains.df[, "maternal:NOD"]
allChains.df$`dam:NZO` <- allChains.df[,"additive:NZO"] + allChains.df[, "maternal:NZO"]
allChains.df$`dam:CAST` <- allChains.df[,"additive:CAST"] + allChains.df[, "maternal:CAST"]
allChains.df$`dam:PWK` <- allChains.df[,"additive:PWK"] + allChains.df[, "maternal:PWK"]
allChains.df$`dam:WSB` <- allChains.df[,"additive:WSB"] + allChains.df[, "maternal:WSB"]
allChains.df$`sire:AJ` <- allChains.df[,"additive:AJ"] - allChains.df[, "maternal:AJ"]
allChains.df$`sire:B6` <- allChains.df[,"additive:B6"] - allChains.df[, "maternal:B6"]
allChains.df$`sire:129` <- allChains.df[,"additive:129"] - allChains.df[, "maternal:129"]
allChains.df$`sire:NOD` <- allChains.df[,"additive:NOD"] - allChains.df[, "maternal:NOD"]
allChains.df$`sire:NZO` <- allChains.df[,"additive:NZO"] - allChains.df[, "maternal:NZO"]
allChains.df$`sire:CAST` <- allChains.df[,"additive:CAST"] - allChains.df[, "maternal:CAST"]
allChains.df$`sire:PWK` <- allChains.df[,"additive:PWK"] - allChains.df[, "maternal:PWK"]
allChains.df$`sire:WSB` <- allChains.df[,"additive:WSB"] - allChains.df[, "maternal:WSB"]
allChains.dam <- allChains.df[, grep(names(allChains.df), pattern="^dam:", value=TRUE)]
allChains.sire <- allChains.df[, grep(names(allChains.df), pattern="^sire:", value=TRUE)]
allChains.fixed <- allChains.df[,c("mu", "inbred", "litternum")]
allChains.add <- as.matrix(allChains.df[, grep(colnames(allChains.df), pattern="^additive[^.]", value=TRUE)])
allChains.mat <- as.matrix(allChains.df[, grep(colnames(allChains.df), pattern="^maternal[^.]", value=TRUE)])
allChains.dam <- as.matrix( allChains.df[, grep(colnames(allChains.df), pattern="^dam[^.]", value=TRUE)])
allChains.sire <- as.matrix(allChains.df[, grep(colnames(allChains.df), pattern="^sire[^.]", value=TRUE)])
allChains.inbred <- as.matrix(allChains.df[, grep(colnames(allChains.df), pattern="^inbred[^.]", value=TRUE)])
allChains.jk <- as.matrix(allChains.df[, grep(colnames(allChains.df), pattern="^v[^.]", value=TRUE)])
allChains.asymm <- as.matrix(allChains.df[, grep(colnames(allChains.df), pattern="^w[^.]", value=TRUE)])
allChains.batch <- as.matrix(allChains.df[, grep(colnames(allChains.df), pattern="^YearMonth[^.]", value=TRUE)])
allChains.batch.1 <- as.matrix(allChains.df[, grep(colnames(allChains.df), pattern="^order[^.]", value=TRUE)])
allChains.sigma <- as.matrix(cbind(allChains.df[, grep(colnames(allChains.df), pattern="^tau-sq-[^.]", value=TRUE)],
`sigma-sq`=allChains.df$`sigma-sq`))
allChains <- as.mcmc(cbind(allChains.fixed, allChains.add, allChains.mat,
allChains.dam, allChains.sire, allChains.inbred,
allChains.jk, allChains.asymm,
allChains.batch, allChains.batch.1,
allChains.sigma))
founder.names <- c("AJ", "B6", "129", "NOD", "NZO", "CAST", "PWK", "WSB")
projection <- as.data.frame(as.matrix(expand.grid(founder.names, founder.names)[,c(2,1)]))
projection$inbred <- ifelse(projection[,1]==projection[,2], 1, 0)
projection.inbred.mat <- diag(projection$inbred)
colnames(projection) <- c("Dam_Founder", "Sire_Founder", "inbred")
projection.matrices <- diallelMatrixMaker(projection, dam.col.name="Dam_Founder", sire.col.name="Sire_Founder")
proj.add <- apply(allChains.add, MARGIN=1, FUN=function(x){projection.matrices$add.mat %*% x})
proj.mat <- apply(allChains.mat, MARGIN=1, FUN=function(x){projection.matrices$mat.mat %*% x})
proj.dam <- apply(allChains.dam, MARGIN=1, FUN=function(x){projection.matrices$dam.mat %*% x})
proj.sire <- apply(allChains.sire, MARGIN=1, FUN=function(x){projection.matrices$sire.mat %*% x})
proj.inbred <- apply(allChains.inbred, MARGIN=1, FUN=function(x){projection.matrices$inbred.mat %*% x})
proj.jk <- apply(allChains.jk, MARGIN=1, FUN=function(x){projection.matrices$jk.mat %*% x})
proj.asymm <- apply(allChains.asymm, MARGIN=1, FUN=function(x){projection.matrices$asymm.mat %*% x})
proj.mu <- matrix(rep(allChains.fixed[,"mu"],each=dim(projection)[1]),nrow=dim(projection)[1])
proj.inbred.overall <- projection.inbred.mat %*% matrix(rep(allChains.fixed[,"inbred"],each=dim(projection)[1]), nrow=dim(projection)[1])
proj.epsilon <- sapply(allChains.sigma[,"sigma-sq"], FUN=function(x){rnorm(n=dim(projection)[1], mean=0, sd=sqrt(x))})
proj.determ.mu <- proj.mu + proj.inbred.overall + proj.add + proj.mat + proj.inbred + proj.jk + proj.asymm
proj.determ <- proj.add + proj.mat + proj.inbred + proj.jk + proj.asymm
# The following projections should be equivalent:
# proj.determ.mu <- proj.mu + proj.inbred.overall + proj.dam + proj.sire + proj.inbred + proj.jk + proj.asymm
# proj.determ <- proj.dam + proj.sire + proj.inbred + proj.jk + proj.asymm
proj.not.add <- proj.determ - proj.add
proj.not.mat <- proj.determ - proj.mat
proj.not.dam <- proj.determ - proj.dam
proj.not.sire <- proj.determ - proj.sire
proj.not.inbred <- proj.determ - proj.inbred
proj.not.inbred.overall <- proj.determ - proj.inbred.overall
proj.not.jk <- proj.determ - proj.jk
proj.not.asymm <- proj.determ - proj.asymm
proj.total <- proj.determ + proj.epsilon
proj.y <- proj.determ.mu + proj.epsilon
ss.add <- NULL
ss.mat <- NULL
ss.dam <- NULL
ss.sire <- NULL
ss.inbred <- NULL
ss.inbred.overall <- NULL
ss.jk <- NULL
ss.asymm <- NULL
ss.epsilon <- NULL
ss.denom <- NULL
for(i in 1:dim(proj.add)[2]){
ss.add[i] <- sum((proj.determ[,i] - proj.not.add[,i])^2)
ss.mat[i] <- sum((proj.determ[,i] - proj.not.mat[,i])^2)
ss.dam[i] <- sum((proj.determ[,i] - proj.not.dam[,i])^2)
ss.sire[i] <- sum((proj.determ[,i] - proj.not.sire[,i])^2)
ss.inbred[i] <- sum((proj.determ[,i] - proj.not.inbred[,i])^2)
ss.inbred.overall[i] <- sum((proj.determ[,i] - proj.not.inbred.overall[,i])^2)
ss.jk[i] <- sum((proj.determ[,i] - proj.not.jk[,i])^2)
ss.asymm[i] <- sum((proj.determ[,i] - proj.not.asymm[,i])^2)
ss.epsilon[i] <- sum((proj.y[,i] - proj.determ.mu[,i])^2)
}
Y.bar <- matrix(rep(apply(proj.determ, MARGIN=2, FUN=sum)/dim(projection)[1], each=dim(projection)[1]), nrow=dim(projection)[1])
ss.denom <- ss.add + ss.mat + ss.inbred + ss.inbred.overall + ss.jk + ss.asymm + ss.epsilon + adjust
additive <- ss.add/ss.denom
mat.dev <- ss.mat/ss.denom
maternal <- ss.dam/ss.denom
paternal <- ss.sire/ss.denom
inbred <- ss.inbred/ss.denom
inbred.overall <- ss.inbred.overall/ss.denom
epistatic.symm <- ss.jk/ss.denom
epistatic.asymm <- ss.asymm/ss.denom
explained <- (ss.add + ss.mat + ss.inbred + ss.jk + ss.asymm)/ss.denom
unexplained <- 1 - explained
varps.all <- cbind(additive, parental.sex=mat.dev, inbred, inbred.overall, symmetric.epistatic=epistatic.symm, asymmetric.epistatic=epistatic.asymm, total.explained=explained, noise=unexplained)
varps.all.table <- cbind(mean=colMeans(varps.all),
lower.bound=HPDinterval(as.mcmc(varps.all))[,1],
upper.bound=HPDinterval(as.mcmc(varps.all))[,2])
varcolors <- c("#A6CEE3", "#B2DF8A", "#FB9A99",
"#FDBF6F", "#CAB2D6", "#D2B48C")
varcolors.y0 <- c(8.5, 7.5, 6.5, 5.5, 4.5, 3.5)
varcolors.y1 <- c(7.5, 6.5, 5.5, 4.5, 3.5, 2.5)
psq <- data.frame(varps.all.table)
par(mar = c(5.1, 13.5, 2.1, 12.1))
rows <- nrow(psq)
plot(y = c(0.5, rows + 0.5), x = c(-0.06, 1.01), yaxs = "i",
xaxs = "i", col = "white", ylab = "", xlab = "Variance Projection",
yaxt = "n")
segments(x0=-0.05, x1=-0.05, y0=varcolors.y0, y1=varcolors.y1, col=varcolors, lwd=15, lend=1)
abline(v = c(0, 0.2, 0.4, 0.6, 0.8, 1), lty = 2, col = "grey")
abline(h = 2.5, lwd = 1.5)
segments(x0 = psq$lower.bound, x1 = psq$upper.bound, y0 = c(rows:1),
y1 = c(rows:1), lwd = 4, lend=1)
axis(side = 2, at = c(rows:1), labels = rownames(psq),
las = 1)
rightlabels <- paste(as.character(sprintf("%.2f", round(100 *
psq$mean, digits = 2))), "% (", as.character(sprintf("%.2f",
round(100 * psq$lower.bound, digits = 2))), "%, ", as.character(sprintf("%.2f",
round(100 * psq$upper.bound, digits = 2))), "%)", sep = "")
axis(side = 4, at = c(rows:1), labels = rightlabels, tick = FALSE,
las = 1)
points(x = psq$mean, y = c(rows:1), pch = 16, col = "black",
cex = 1.5)
points(x = psq$mean, y = c(rows:1), pch = 16, col = "white",
cex = 0.5)
Table of results
print(xtable(psq, digits=4), type="html")
|
|
mean
|
lower.bound
|
upper.bound
|
|
additive
|
0.0914
|
0.0558
|
0.1331
|
|
parental.sex
|
0.0569
|
0.0286
|
0.0838
|
|
inbred
|
0.0035
|
0.0001
|
0.0071
|
|
inbred.overall
|
0.0142
|
0.0087
|
0.0209
|
|
symmetric.epistatic
|
0.0054
|
0.0006
|
0.0147
|
|
asymmetric.epistatic
|
0.0183
|
0.0051
|
0.0325
|
|
total.explained
|
0.1755
|
0.1190
|
0.2378
|
|
noise
|
0.8245
|
0.7622
|
0.8810
|
Plot dam_strain and sire_strain VarPs
varps.all.damsire <- cbind(dam.strain=maternal, sire.strain=paternal, inbred, inbred.overall, symmetric.epistatic=epistatic.symm, asymmetric.epistatic=epistatic.asymm, total.explained=explained, noise=unexplained)
varps.all.damsire.table <- cbind(mean=colMeans(varps.all.damsire),
lower.bound=HPDinterval(as.mcmc(varps.all.damsire))[,1],
upper.bound=HPDinterval(as.mcmc(varps.all.damsire))[,2])
## MCMCglmm model variance component estimates
varcolors <- c("#781112", "#131E3A", "#FB9A99",
"#FDBF6F", "#CAB2D6", "#D2B48C")
varcolors.y0 <- c(8.5, 7.5, 6.5, 5.5, 4.5, 3.5)
varcolors.y1 <- c(7.5, 6.5, 5.5, 4.5, 3.5, 2.5)
psq <- data.frame(varps.all.damsire.table)
par(mar = c(5.1, 13.5, 2.1, 12.1))
rows <- nrow(psq)
plot(y = c(0.5, rows + 0.5), x = c(-0.06, 1.01), yaxs = "i",
xaxs = "i", col = "white", ylab = "", xlab = "Variance Projection",
yaxt = "n")
segments(x0=-0.05, x1=-0.05, y0=varcolors.y0, y1=varcolors.y1, col=varcolors, lwd=15, lend=1)
abline(v = c(0, 0.2, 0.4, 0.6, 0.8, 1), lty = 2, col = "grey")
abline(h = 2.5, lwd = 1.5)
segments(x0 = psq$lower.bound, x1 = psq$upper.bound, y0 = c(rows:1),
y1 = c(rows:1), lwd = 4, lend=1)
axis(side = 2, at = c(rows:1), labels = rownames(psq),
las = 1)
rightlabels <- paste(as.character(sprintf("%.2f", round(100 *
psq$mean, digits = 2))), "% (", as.character(sprintf("%.2f",
round(100 * psq$lower.bound, digits = 2))), "%, ", as.character(sprintf("%.2f",
round(100 * psq$upper.bound, digits = 2))), "%)", sep = "")
axis(side = 4, at = c(rows:1), labels = rightlabels, tick = FALSE,
las = 1)
points(x = psq$mean, y = c(rows:1), pch = 16, col = "black",
cex = 1.5)
points(x = psq$mean, y = c(rows:1), pch = 16, col = "white",
cex = 0.5)
Table of results
print(xtable(psq, digits=4), type="html")
|
|
mean
|
lower.bound
|
upper.bound
|
|
dam.strain
|
0.1304
|
0.0849
|
0.1843
|
|
sire.strain
|
0.0180
|
0.0044
|
0.0335
|
|
inbred
|
0.0035
|
0.0001
|
0.0071
|
|
inbred.overall
|
0.0142
|
0.0087
|
0.0209
|
|
symmetric.epistatic
|
0.0054
|
0.0006
|
0.0147
|
|
asymmetric.epistatic
|
0.0183
|
0.0051
|
0.0325
|
|
total.explained
|
0.1755
|
0.1190
|
0.2378
|
|
noise
|
0.8245
|
0.7622
|
0.8810
|