Skip to contents

Find the AIC and BIC for a model fitted using SpATS

Source: R/icREML.R
ic_reml_spt.Rd

This function calculates the AIC and BIC for a model fitted in SpATS following the methodology proposed by Verbyla (2019).

Usage

ic_reml_spt(model, scale = 1, k = 2, label = "spats")

Arguments

model

A model fitted using SpATS.

scale

A scalar to scale the variance matrix of the estimated fixed effects (to ensure numerical stability of a log-determinant). Default value is 1.

k

An integer value to round ratios when identifying boundary variance parameters. Default value is 2.

label

A string to label the model. Default value is "spats".

Value

A data frame. The data frame has the following components

  • model : the name of the models

  • loglik : the full log-likelihood for each model

  • p : the number of fixed effects parameters for each model

  • q : the number of (non-zero) variance parameters for each model.

  • b : the number of variance parameters that are fixed or on the boundary. These parameters are not counted in the AIC or BIC.

  • AIC : the AIC for each model

  • BIC : the BIC for each model

  • logdet : the log-determinant used in adjusting the residual log-likelihood for each model

References

Verbyla, A. P. (2019). A note on model selection using information criteria for general linear models estimated using REML. Australian & New Zealand Journal of Statistics, 61(1), 39-50.

Author

Johan Aparicio

Examples

# \donttest{
library(SpATS)
library(agriutilities)
data(wheatdata)
wheatdata$R <- as.factor(wheatdata$row)
wheatdata$C <- as.factor(wheatdata$col)

m1 <- SpATS(
  response = "yield",
  spatial = ~ PSANOVA(col, row, nseg = c(10, 20), nest.div = 2),
  genotype = "geno",
  genotype.as.random = TRUE,
  fixed = ~ colcode + rowcode,
  random = ~ R + C,
  data = wheatdata,
  control = list(tolerance = 1e-03, monitoring = 0)
)

m2 <- SpATS(
  response = "yield",
  spatial = ~ PSANOVA(col, row, nseg = c(10, 20), nest.div = 2),
  genotype = "geno",
  genotype.as.random = TRUE,
  fixed = ~colcode,
  random = ~ R + C,
  data = wheatdata,
  control = list(tolerance = 1e-03, monitoring = 0)
)

rbind.data.frame(
  ic_reml_spt(m1, label = "colcode_rowcode"),
  ic_reml_spt(m2, label = "colcode_no_rowcode")
)
#>                model res_loglik full_loglik p q b      AIC      BIC   logdet
#> 1    colcode_rowcode  -1520.643   -1551.460 8 8 1 3134.920 3195.705 61.63441
#> 2 colcode_no_rowcode  -1527.611   -1556.439 7 8 1 3142.878 3199.864 57.65614

rbind.data.frame(
  h_cullis_spt(m1),
  h_cullis_spt(m2)
)
#>   trait  H2Cullis H2Oakey reBLUP_avg vdBLUP_avg  PEV_avg    var_G
#> 1 yield 0.7726571    0.77  0.7654360   1192.740 615.3125 2623.218
#> 2 yield 0.7691996    0.77  0.7620109   1207.598 622.6054 2616.109
# }