proLik | R Documentation |
Estimation, checking, and plotting of profile likelihoods and objects of
class proLik
from a mixed model in ASReml-R.
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 class 'proLik' is(x) ## S3 method for class 'proLik' plot( x, CL = TRUE, alpha = NULL, type = "l", main = NULL, xlab = NULL, ylab = NULL, ... )
full.model |
An |
component |
A character (alternatively for |
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 |
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 |
CL |
A logical indicating whether a lines at the confidence limits are to be drawn when plotting. |
type, main, xlab, ylab |
See arguments to |
... |
other arguments to |
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.
An S3 object of class “proLik” containing:
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
value along the sampling interval for which the
component
was constrained
approximate Upper Confidence Limit
approximate Lower Confidence Limit
the component for which the profile likelihood surface has been constructed
the alpha level at which the confidence limits have been calculated
May be unfeasible to estimate profile likelihoods for complex models with many variance components
aiFun
## Not run: 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) ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.