nadiv (version 2.11)

proLik: Estimates the profile likelihood of a random effect

Description

When a mixed model is run in ASReml-R, this function can estimate the profile likelihood of variance components from the output.

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

Arguments

full.model

An asreml model object

component

A character indicating for which variance component the profile likelihood will be constructed. Must be an object in full.model$gammas.

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 sample values and interpolated values, for determining the upper and lower limits of the Confidence Interval. Actual sample points will be no more than this distance from the true value of the estimate.

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.

Value

• 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

Warning

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

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.

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

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 ~ 1, 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.proLik(profileA) 
      plot.proLik(profileD)
      plot.proLik(profileE)