nadiv (version 2.17.2)

proLik: Profile Likelihoods

Description

Estimation, checking, and plotting of profile likelihoods and objects of class proLik from a mixed model in ASReml-R.

Usage

proLik(
  full.model,
  component,
  G = TRUE,
  negative = FALSE,
  nsample.units = 3,
  nse = 3,
  alpha = 0.05,
  tolerance = 0.001,
  parallel = FALSE,
  ncores = getOption("mc.cores", 2L)
)
proLik4(
  full.model,
  component,
  G = TRUE,
  negative = FALSE,
  nsample.units = 3,
  nse = 3,
  alpha = 0.05,
  tolerance = 0.001,
  parallel = FALSE,
  ncores = getOption("mc.cores", 2L)
)
# S3 method for proLik
is(x)
# S3 method for proLik
plot(
  x,
  CL = TRUE,
  alpha = NULL,
  type = "l",
  main = NULL,
  xlab = NULL,
  ylab = NULL,
  ...
)

Value

An S3 object of class “proLik” containing:

lambdas

negative log-Likelihood ratio test statistic. Estimated from the log-Likelihood of the full.model and the log-Likelihood of the model with the component constrained to a value in the sampling interval

var.estimates

value along the sampling interval for which the component was constrained

UCL

approximate Upper Confidence Limit

LCL

approximate Lower Confidence Limit

component

the component for which the profile likelihood surface has been constructed

alpha

the alpha level at which the confidence limits have been calculated

Arguments

full.model

An asreml model object

component

A character (alternatively for proLik4 this could also be a formula) indicating for which variance component the profile likelihood will be constructed. For proLik, must be an object in full.model$gammas. For proLik4, must be an object in full.model$vparameters. To specify as a formula, components are written as Vx, where “x” is a number between 1 and length(full.model$vparameters) (e.g., component = ~ V1).

G

Logical indicating whether component is part of the G structure. If the component is part of the R structure, G = FALSE.

negative

Logical indicating whether or not the component can take on a negative value (i.e., a covariance)

nsample.units

Number of sample units to be used in constructing the area between the confidence limits for the profile likelihood

nse

Number of standard errors on either side of the estimate, over which the confidence limits should be evaluated.

alpha

The critical value for determining the Confidence Interval

tolerance

Acceptable distance between actual and desired minimization criteria for determining the upper and lower limits of the confidence interval. See Details for more

parallel

A logical indicating whether or not parallel processing will be used. Note, may only be available for Mac and Linux operating systems.

ncores

Argument indicating number of cpu units to use. Default is all available.

x

Object of class proLik.

CL

A logical indicating whether a lines at the confidence limits are to be drawn when plotting.

type, main, xlab, ylab

See arguments to plot.

...

other arguments to plot.

Warning

May be unfeasible to estimate profile likelihoods for complex models with many variance components

Author

Details

For the negative argument, this should be used if the profile likelihood of a covariance component is to be constructed.

parallel = TRUE should only be used on Linux or Mac OSes (i.e., not Windows).

The function uses the optimize function to obtain the approximate confidence limits. Therefore, nse should be carefully thought about beforehand when running the function. Increasing this value will ensure the confidence limits are contained by the search space, but at a cost to time.

tolerance is expressed in terms of the alpha value desired. Convergence to a confidence limit will only be achieved when the constrained value proposed for the confidence limit yields a likelihood ratio test statistic that is no worse than alpha + tolerance.

If negative is FALSE, and the lower bound of the sampling interval extends beyond zero, this will instead be set to effectively zero.

Obtaining the profile likelihood for residual variances may necessitate explicitly specifying the R structure of the ASReml model. See example below.

See Also

aiFun

Examples

  if (FALSE) {
    library(asreml)
    ginvA <- asreml.Ainverse(warcolak[, c(1,3,2)])$ginv
    ginvD <- makeD(warcolak[,1:3])$listDinv
    warcolak$IDD <- warcolak$ID
    warcolak.mod <- asreml(trait1 ~ sex, random = ~ ped(ID) + giv(IDD), 
	ginverse = list(ID = ginvA, IDD = ginvD), rcov = ~ idv(units), data = warcolak) 
    summary(warcolak.mod)$varcomp
    profileA <- proLik(full.model = warcolak.mod, component = "ped(ID)!ped", 
        G = TRUE, negative = FALSE, nsample.units = 3, nse = 3)
    profileA
    profileD <- proLik(warcolak.mod, component = "giv(IDD).giv", 
	G = TRUE, negative = FALSE, nsample.units = 3, nse = 3)
    profileE <- proLik(warcolak.mod, component = "R!units.var", G = FALSE, negative = FALSE)

    x11(w = 6, h = 8)
    par(mfrow = c(3,1))
      plot(profileA) 
      plot(profileD)
      plot(profileE)
   }