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

covaruber / sommer