nadiv (version 2.9)

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 = 4, nse = 2, alpha = 0.05, threshold = 0.001,
       parallel = FALSE, ncores = getOption("cores"))

Arguments

full.model

An asreml model object

component

A character vector 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 profile likelihood

nse

Number of standard errors on either side of the estimate, over which the profile likelihood should be constructed

alpha

The critical value for determining the Confidence Interval

threshold

Acceptable distance, between actual sample values and interpolated values, for determining the upper and lower limits of the Confidence Interval. Expressed as a percentage of the parameter estimate. Actual sample points will be no more than this distanc

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.

If parallel = TRUE then the pacakge multicore must be loaded.

The function uses a grid search to obtain the approximate confidence limits. Therefore, nsample.units should be carefully thought about beforehand when running the function. Increasing this value will ensure a smoother surface, but at cost to time. The function weights this quantity to ensure more samples for regions of the profile likelihood greater than the estimate. The lower bound of the sampling interval to the estimate is made up of 2 times nsample.units, whereas the sampling interval from the estimate to the upper bound is 3 times nsample.units.

Similarly for nse, the function will include twice as many standard errors from the estimate to the upper bound than it does from the estimate to the lower bound. If negative is FALSE, and the lower bound of the sampling interval extends beyond zero, this will instead be set to 0.

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), 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(full.model = warcolak.mod, component = "giv(IDD).giv", 
	G = TRUE, negative = FALSE, nsample.units = 3, nse = 3)

    x11(w = 6, h = 8)
    par(mfrow = c(2,1))
      plot.proLik(profileA) 
      plot.proLik(profileD)