ic_reml_spt: Find the AIC and BIC for a model fitted using SpATS

In AparicioJohan/agriutilities: Utilities for Data Analysis in Agriculture

Find the AIC and BIC for a model fitted using SpATS

Description

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

Author(s)

Johan Aparicio

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.

Examples

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")
)

rbind.data.frame(
  h_cullis_spt(m1),
  h_cullis_spt(m2)
)

AparicioJohan/agriutilities documentation built on Jan. 20, 2025, 7:53 a.m.