mean comparisons with spatial model

[priviere]

is mean comparisons method ok ?

In order to assess mean comparisons, I used the agricolae::LSD.test function

When doing ?LSD.test, in example, we can see

# Old version LSD.test()
df<-df.residual(model)
MSerror<-deviance(model)/df
out <- with(sweetpotato,LSD.test(yield,virus,df,MSerror))

So I decided to follow the same logic in model_spatial:

s = summary.SpATS(m, which = "dimensions")
deviance_model = deviance(m)
df_residual = as.numeric(s$p.table.dim["Residual", "Effective"])
MSerror = deviance_model / df_residual

The degree of freedom of the residuals are the effective degree of freedom

And then in mean_comparisons.check_model_spatial:

  lsd = agricolae::LSD.test(
    y = data[,variable],
    trt = data$germplasm,
    DFerror = x$spatial$model$df_residual, 
    MSerror = x$spatial$model$MSerror, 
    alpha = alpha, 
    p.adj = p.adj
    )

I hope it is correct ...

[priviere]

Look at confidence interval for BLUPs
https://stackoverflow.com/questions/13847936/in-r-plotting-random-effects-from-lmer-lme4-package-using-qqmath-or-dotplot

[priviere]

It may be different if genotype are random or not

• if random: see with BLUP
• if fixed: the first post with LSD may be ok ?

[priviere]

vcov in SpATS output may be useful as, according to help, it is nverse of the coefficient matrix of the mixed models equations. The inverse is needed for the computation of standard errors. For computational issues, the inverse is returned as a list: C11_inv corresponds to the coefficient matrix associated with the genotype; C22_inv corresponds to the coefficient matrix associated with the spatial, the fixed and the random components; and C12_inv and C21_inv correspond to the “combination” of both.

C11_inv can be used to get confidence interval in order then to do mean comparisons

[priviere]

library(SpATS)
data(wheatdata)
summary(wheatdata)

# Create factor variable for row and columns
wheatdata$R <- as.factor(wheatdata$row)
wheatdata$C <- as.factor(wheatdata$col)

m0 <- SpATS(response = "yield", spatial = ~ SAP(col, row, nseg = c(10,20), degree = 3, pord = 2), 
            genotype = "geno", fixed = ~ colcode + rowcode, random = ~ R + C, data = wheatdata, 
            control =  list(tolerance = 1e-03))

pre = predict.SpATS(m0, which = "geno")
pre$geno <- factor(pre$geno, levels = pre$geno[order(pre$predicted.values)])
pre$lower = pre$predicted.values-pre$standard.errors
pre$upper = pre$predicted.values+pre$standard.errors
ggplot(pre, aes(y = predicted.values, x = geno)) + geom_point() + geom_errorbar(aes(ymin = lower, ymax = upper)) + coord_flip()

[priviere]

Regarding LSD.test, it gives coherent results with the following code:

s = summary.SpATS(m0, which = "dimensions")
deviance_model = deviance(m0)
df_residual = as.numeric(s$p.table.dim["Residual", "Effective"])
MSerror = deviance_model / df_residual

lsd = agricolae::LSD.test(
  y = pre$predicted.values,
  trt = pre$geno,
  DFerror = df_residual, 
  MSerror = MSerror, 
  alpha = 0.05,
  p.adj = "none"
)

Output from BLUPs ggplot and groups with LSD.test are coherent !
Nevertheless, this seems strange to do that. When doing test, it appears the following bug from agricolae::LSD.test(): Error in if (lastC(M[i]) != letras[k]) M[i] <- paste(M[i], letras[k], : missing value where TRUE/FALSE needed
It is because there are too many groups and the function agricolae::orderPvalue() bugs because letras object, which is a vector containing the letter, is too short.
On trick to handle this is to settle alpha = 1 in mean_comparisons(): out_mean_comparisons_spatial = mean_comparisons(out_check_spatial, alpha = 1, p.adj = "bonferroni")

[priviere]

code updated in v0.25