onset.test {OnAge}

Tests the null hypothesis that the age at the onset of senescence is the same in two groups

Description

The function takes as input a log-likelihood function defining a statistical model of the senescence and two datasets. It performs a likelihood ratio test of the null hypothesis that the age at the onset of senescence is the same in the two groups.

Usage

onset.test(ll, data1, data2, search.range, CI.lvl=0.95, 
    tol=.Machine$double.eps^0.25, warn=FALSE, do.plot=FALSE, 
    plot.file=NULL, grid.len=100)

Arguments

ll

A function computing the log-likelihood of data under a statistical model of the senescence.

data1

A data.frame object containing the data for the first group.

data2

A data.frame object containing the data for the second group.

search.range

The range over which the log-likelihood should be maximized with respect to the age at the onset of senescence.

CI.lvl

An optional confidence level to obtain a confidence interval for each onset estimate. No confidence interval is computed if NA. Default 0.95.

tol

A tolerance, to be passed to the optimize function. Optional.

warn

A logical variable indicating whether the function should issue a warning when the computed likelihood ratio is negative (which may happen for numerical reasons). Optional.

do.plot

A logical variable indicating whether the function should produce a plot of the log-likelihood against age at the onset of senescence (with vertical dotted lines indicating the confidence interval if CI.lvl is not NA). Optional (default FALSE).

plot.file

If not NULL and do.plot is TRUE, the function makes the plot in a pdf file with name plot.file. Optional (default NULL).

grid.len

Integer value. If do.plot is TRUE, specifies the number of onset values for which the log-likelihood is evaluated between search.range[1] and search.range[2]. Optional (default 100).

Value

A list with elements:

pv

A p-value for the tested null hypothesis.

est.1

A maximum likelihood estimate for the age at the onset of senescence in the first group.

est.2

A maximum likelihood estimate for the age at the onset of senescence in the second group.

est.joint

A maximum likelihood estimate for the age at the onset of senescence in the merged groups.

CI.1

A confidence interval for the age at the onset of senescence in the first group, with confidence level CI.lvl (NA if CI.lvl is NA).

CI.2

A confidence interval for the age at the onset of senescence in the second group, with confidence level CI.lvl (NA if CI.lvl is NA).

joint.CI

A confidence interval for the age at the onset of senescence in the total group, with confidence level CI.lvl (NA if CI.lvl is NA).

lh0

The log-likelihood maximized under the null hypothesis.

lh1

The log-likelihood maximized under the alternative hypothesis.

llr

The likelihood ratio statistic.

cvg.ok

A logical variable indicating whether the computed likelihood ratio was negative (the returned value is 0 in this case).

Examples

if(requireNamespace("lme4", quietly=TRUE)) {
    data(RoeDeerMassData)
    RoeDeerMassData$ID <- factor(RoeDeerMassData$ID)
    RoeDeerMassData$cohort <- factor(RoeDeerMassData$cohort)

    dataFCH <- RoeDeerMassData[RoeDeerMassData$sex%in%"F"&
                                 RoeDeerMassData$population%in%"CH", ]
    dataMCH <- RoeDeerMassData[RoeDeerMassData$sex%in%"M"&
                                 RoeDeerMassData$population%in%"CH", ]
    dataFTF <- RoeDeerMassData[RoeDeerMassData$sex%in%"F"&
                                 RoeDeerMassData$population%in%"TF", ]
    dataMTF <- RoeDeerMassData[RoeDeerMassData$sex%in%"M"&
                                 RoeDeerMassData$population%in%"TF", ]

    ## b1: function for piecewise regression (transforms x into 0 before bp)
    b1 <- function(x, bp) ifelse(x < bp, 0, x - bp)

    ## Use this function to define the model in which the differential
    ## onset hypothesis is tested.
    ll <- function(thr, dataIn){
        logLik(lme4::lmer(body.mass ~ b1(age, thr) + age.at.last.capture + 
            last.year.of.capture + (1|ID) + (1|cohort), data=dataIn, 
            REML="FALSE"))
    }

    search.range <- c(6,12) # data not available before 6 years old
    search.range.TF <- search.range.CH <- search.range

    ## Not run: 
      # if you don't run the example in RStudio, create a new graphic 
      # window with suitable dimensions
      if(.Platform$GUI!="RStudio") {dev.new(width=15, height=5)}
      # all graphs are square
      par(pty="s")
      
      # test for Trois Fontaines
      res.tf <- onset.test(ll, dataFTF, dataMTF, search.range.TF, do.plot=TRUE)
      
      if(.Platform$GUI!="RStudio") {dev.new(width=15, height=5)}
      par(pty="s")
      
      # test for Chizé
      res.ch <- onset.test(ll, dataFCH, dataMCH, search.range.CH, do.plot=TRUE)
      
      cat(sprintf("p-value for differential age at onset is %g in 
          Trois Fontaines, %g in Chizé", res.tf$pv, res.ch$pv))
    
## End(Not run)
    
    # a shorter executable example for Trois Fontaines
    set.seed(40)
    shortdataFTF <- dataFTF[sample(1:NROW(dataFTF), 200), ]
    shortdataMTF <- dataMTF[sample(1:NROW(dataMTF), 100), ]
    res.tf.short <- onset.test(ll, shortdataFTF, shortdataMTF, search.range.TF, 
      do.plot=FALSE, CI.lvl=NA)
  }