R: Estimate confidence intervals using profile likelihood
Profile.like.CI
R Documentation
Estimate confidence intervals using profile likelihood
Description
profile likelihood is used to estimate 95
Usage
Profile.like.CI(DIST, TIME, GRAD, meserr1 = 0, meserr2 = 0, like, par,
MODEL, test.values.par1, test.values.par2, p_starting="NULL")
Arguments
DIST
vector of Euclidean distances for sister pair dataset
TIME
vector of evolutionary ages (i.e. node ages ) for sister pair dataset
GRAD
vector of gradient values (i.e. any continuous variable) for sister pair dataset (see Details)
meserr1
a list of measurement errors that correspond to the first of each species in a sister pair. Order of sister pairs is the same as for DIST.
meserr2
a list of measurement errors that correspond to the second of each species in a sister pair. Order of sister pairs is the same as for DIST.
like
loglike at the MLE as returned by model.test.sisters
par
a vector containing the maximum likelihood estimates of model parameters. These should be in the order indicated in sisterContinuous.
MODEL
The name of the gradient model to perform profile likelihood on. Currently implemented models are "BM_linear", "OU_linear_beta", and "OU_linear"
test.values.par1
a vector of values to calculate likelihoods for parameter 1
test.values.par2
a vector of values to calculate likelihoods for parameter 2
p_starting
List of starting values for the model. If starting=list("NULL"), the built-in starting parameters are used, and is generally recommended.
Details
This function uses profile likelihood to estimate confidence for select parameters. Likelihood surfaces often have ridges (e.g. the "OU_linear" model), and the resulting confidence intervals are not always symmetric around the MLE.
Profile likelihood generates confidence intervals appropriate in such cases. However, the code is computationally demanding.
Currently, profile likelihood is only implemented for the two parameters of BM_linear and for the slope parameters of OU_linear_beta (1 parameter) and OU_linear (2 parameters).
Value
Returns a list with the following elements:
profile.likelihoods_par1, 2, 3 etc For each parameter, a matrix showing the range of values tested (test.value), the log likelihoods of each value in the range (logLike), the difference in likelihood from the MLE and each value (logLikeDifference). The final column (CI_range) gives a 1 if the value was less than 1.92 loglikelihood units below the MLE, and thus outside the 95
model The model used
MLE_par1 The MLE for parameter 1
CI_par1 The lower and upper 95
warnings_par1 A warning is returned only if the lower or upper limit of the CI has not been reached by the range of tested values. Otherwise, returns NA
Author(s)
Jason T. Weir
See Also
bootstrap.test
Examples
## Not run:
###This example uses the syllable dataset for oscine songbirds Weir & Wheatcroft 2011
data(bird.syllables)
attach(bird.syllables)
#STEP 1 Correct Euclidean distances for sampling and measurement bias
DIST_cor <- MScorrection(nA=bird.syllables$number_individuals_Species1,
nB=bird.syllables$number_individuals_Species2,
VarA=bird.syllables$Species1_PC2_var,
VarB=bird.syllables$Species2_PC2_var,
DIST_actual=abs(bird.syllables$Species1_PC2_mean -
bird.syllables$Species2_PC2_mean))
#STEP 2 Test all models on oscines only (in which song has a strong
#culturally transmitted component)
DIST <- subset(DIST_cor, subset = (bird.syllables$Suboscine == "oscine"))
TIME <- subset(bird.syllables$TIME,subset = (bird.syllables$Suboscine == "oscine"))
GRAD <- subset(bird.syllables$GRAD,
subset = (bird.syllables$Suboscine == "oscine"))
FIT5 <- model.test.sisters(DIST=DIST, TIME=TIME, GRAD=GRAD, models=models)
#The best fit model in FIT5 is BM_linear in which tropical species have a
#much slower rate than temperate species.
#STEP 3 run the profile likelihood
Profile.like.CI(DIST=DIST, TIME=TIME, GRAD=GRAD, meserr1 = 0, meserr2 = 0,
like=FIT5[1,2], par=c(FIT5[5,2], FIT5[6,2]), MODEL="BM_linear", MULT=1,
test.values.par1 = c((0:100)*0.001), test.values.par2 = c((33:100)*0.0001),
p_starting="NULL")
## End(Not run)#end dontrun