sommer
sommer copied to clipboard
Published
20 hours ago •
covaruber
covaruber
Different varcomp estimates between sommer and ASReml-R when it is bounded
Hi Giovanny,
I found different results between sommer and ASReml-R for a simple alpha-design analysis. The mod$beta already includes intercept and the predict.mmer function results in the same values as the mod$beta. Also, the varcomp estimates are different from ASReml-R. Perhaps, this is the reason that the results are different.
I provide the code (and results) below.
Upload the package
library(sommer)
library(data.table)
library(agridat)
library(asreml)
options("scipen" = 100,"digits" = 4 ) # set numbering format
The data
df <- buntaran.wheat
Env <- levels(factor(paste(df$zone, df$loc, sep = "_")))
i = Env[1]
i
[1] "middle_07bm27"
Analysis with sommer
Edat <- droplevels(subset(df, Env == i))
mod.1 <- mmer(fixed = yield ~ gen-1,
random = ~ rep + rep:alpha,
rcov = ~ units,
data = Edat)
summary(mod.1)$varcomp
VarComp VarCompSE Zratio Constraint
rep.yield-yield -710 4346 -0.1634 Positive
rep:alpha.yield-yield 0 7804 0.0000 Positive
units.yield-yield 49331 13338 3.6985 Positive
mod.1$Beta
Trait Effect Estimate
1 yield genG22455 783.7
2 yield genG23286 839.2
3 yield genG24054 1003.2
4 yield genG24521 801.4
5 yield genG25512 896.3
6 yield genG25965 852.1
7 yield genG26742 897.4
8 yield genG26777 829.7
9 yield genG27125 934.2
10 yield genG27130 903.1
11 yield genG27543 425.4
12 yield genG27546 800.7
13 yield genG27548 798.6
14 yield genG27590 986.2
15 yield genG27592 885.2
16 yield genG27593 872.0
17 yield genG27600 910.9
18 yield genG27605 832.1
19 yield genG27669 888.3
20 yield genG28128 805.4
21 yield genG28209 911.1
22 yield genG28949 970.0
23 yield genG28950 841.8
Predict BLUE - sommer
blue <- predict.mmer(mod.1, classify = "gen") # get the lsmeans and vcov
blue.1 <- data.table(blue$pvals)[, c(2:4)] # grab the lsmeans
names(blue.1) <- c("gen", "Yield_adj", "se") # rename columns
blue.1[ , ':='(smith.w = diag(solve(mod.1$VarBeta)))] # calculate the Smith's weight
blue.1
gen Yield_adj se smith.w
1: G22455 783.7 154.5 0.00004260
2: G23286 839.2 108.9 0.00008622
3: G24054 1003.2 154.5 0.00004260
4: G24521 801.4 107.4 0.00008930
5: G25512 896.3 154.5 0.00004260
6: G25965 852.1 108.9 0.00008622
7: G26742 897.4 154.5 0.00004260
8: G26777 829.7 108.9 0.00008622
9: G27125 934.2 220.3 0.00002078
10: G27130 903.1 108.9 0.00008622
11: G27543 425.3 220.3 0.00002078
12: G27546 800.8 154.5 0.00004260
13: G27548 798.5 126.6 0.00006338
14: G27590 986.2 220.3 0.00002078
15: G27592 885.2 154.5 0.00004260
16: G27593 872.0 108.9 0.00008622
17: G27600 910.9 220.3 0.00002078
18: G27605 832.1 108.9 0.00008622
19: G27669 888.3 154.5 0.00004260
20: G28128 805.4 109.4 0.00008519
21: G28209 911.1 220.3 0.00002078
22: G28949 970.0 154.5 0.00004260
23: G28950 841.8 108.9 0.00008622
gen Yield_adj se smith.w
Analysis with ASReml-R
mod.1.asr <- asreml(fixed = yield ~ gen,
random = ~ rep + rep:alpha,
data = Edat)
mod.1.asr <- update.asreml(mod.1.asr)
summary(mod.1.asr)$varcomp
component std.error z.ratio bound %ch
rep 0.004965 NA NA B 0
rep:alpha 0.078505 NA NA B 0
units!R 49065.721485 11408 4.301 P 0
sm <- summary(mod.1.asr, coef =T)
sm$coef.fixed
solution std error z.ratio
gen_G22455 0.00 NA NA
gen_G23286 69.61 191.8 0.36288
gen_G24054 219.56 221.5 0.99121
gen_G24521 36.51 191.8 0.19033
gen_G25512 112.62 221.5 0.50843
gen_G25965 82.49 191.8 0.43003
gen_G26742 113.71 221.5 0.51335
gen_G26777 60.13 191.8 0.31345
gen_G27125 150.50 271.3 0.55474
gen_G27130 133.57 191.8 0.69631
gen_G27543 -339.50 271.3 -1.25144
gen_G27546 17.09 221.5 0.07714
gen_G27548 27.42 202.2 0.13560
gen_G27590 202.55 271.3 0.74663
gen_G27592 101.57 221.5 0.45853
gen_G27593 102.40 191.8 0.53381
gen_G27600 127.25 271.3 0.46904
gen_G27605 62.56 191.8 0.32612
gen_G27669 104.61 221.5 0.47228
gen_G28128 31.11 191.8 0.16218
gen_G28209 127.42 271.3 0.46968
gen_G28949 186.36 221.5 0.84133
gen_G28950 72.25 191.8 0.37665
(Intercept) 774.26 156.6 4.94326
Predict BLUE - ASReml-R
blue.asr <- predict(mod.1.asr, classify = "gen", vcov = TRUE) # get the lsmeans and vcov
blue.1.asr <- data.table(blue.asr$pvals)[, c(1:3)] # grab the lsmeans
names(blue.1.asr) <- c("gen", "Yield_adj", "se") # rename columns
blue.1.asr[ , ':='(smith.w = diag(solve(blue.asr$vcov)))] # calculate the Smith's weight
blue.1.asr
gen Yield_adj se smith.w
1: G22455 774.3 156.6 0.00004076
2: G23286 843.9 110.8 0.00008152
3: G24054 993.8 156.6 0.00004076
4: G24521 810.8 110.8 0.00008152
5: G25512 886.9 156.6 0.00004076
6: G25965 856.8 110.8 0.00008152
7: G26742 888.0 156.6 0.00004076
8: G26777 834.4 110.8 0.00008152
9: G27125 924.8 221.5 0.00002038
10: G27130 907.8 110.8 0.00008152
11: G27543 434.8 221.5 0.00002038
12: G27546 791.3 156.6 0.00004076
13: G27548 801.7 127.9 0.00006114
14: G27590 976.8 221.5 0.00002038
15: G27592 875.8 156.6 0.00004076
16: G27593 876.7 110.8 0.00008152
17: G27600 901.5 221.5 0.00002038
18: G27605 836.8 110.8 0.00008152
19: G27669 878.9 156.6 0.00004076
20: G28128 805.4 110.8 0.00008152
21: G28209 901.7 221.5 0.00002038
22: G28949 960.6 156.6 0.00004076
23: G28950 846.5 110.8 0.00008152
gen Yield_adj se smith.w
Best wishes, Harimurti
I still haven't fixed this issue with mmer but just for the record so users are aware that this bug is not present if using the mmec() function which returns the same estimates. The two stage approach can be done successfully with mmec.
Cheers, Eduardo
