Fit a COD GMRF model using weighted least squares

Usage

codls(tr1, logtau = NULL, profcontrol = list(), weights = NULL, ncpu = 1)

Arguments

tr1: Phylogenetic tree in ape::phylo format
logtau: Precision parameter. If NULL, will invoke `tauprofile` to find best value.
profcontrol: Optional list of arguments passed to `tauprofile`
weights: An optional vector (named or unnamed) of sample weights for each tip in the input tree
ncpu: Integer number of cpu's to use if parallel calculation of tau profile is desired

Value

A COD GMRF model fit

Examples

# A simple example that does not have population structure 
set.seed( 1111 )
tr <- ape::rtree( 100 )
f <- codls(tr)
#>    logtau     loss optimal
#> 1    -4.0 839.6102        
#> 2     0.1 372.3983        
#> 3     4.2 345.3068        
#> 4     8.3 345.2809        
#> 5    12.4 345.2809        
#> 6    16.5 345.2809     ***
#> 7    20.6      Inf        
#> 8    24.7      Inf        
#> 9    28.8      Inf        
#> 10   32.9      Inf        
#> 11   37.0      Inf        
coef(f) |> head() 
#> [1] -6.139952e-14 -5.284621e-14 -5.814045e-14 -5.348977e-14 -8.355911e-14
#> [6] -8.360543e-14
summary(f)
#>  Genealogical placement GMRF model fit 
#> 
#> Phylogenetic tree with 100 tips and 99 internal nodes.
#> 
#> Tip labels:
#>   t13, t36, t48, t79, t65, t66, ...
#> 
#> Rooted; includes branch length(s).
#> Range of coefficients: 
#> -1.21429158048594e-13 1.49802723784303e-13
#> Precision parameter (log tau): 16.5 
#> 
#>   logprecision       RMSCLO Neff
#> 1         16.5 7.152838e-14  100
#> 
if (FALSE) { # \dontrun{
plot(f)
} # }

# This example has population structure 
tr0 = rcoal(20); tr0$edge.length <- .01*tr0$edge.length 
tr1 = rcoal(80); 
dx <- (max(node.depth.edgelength( tr1 ))-max(node.depth.edgelength( tr0 )))
tr0$root.edge <- dx
tr <- bind.tree(tr0,tr1, position = dx)
f <- codls(tr)
#>    logtau     loss optimal
#> 1    -4.0 3743.248        
#> 2     0.1 3210.785     ***
#> 3     4.2 3228.344        
#> 4     8.3 3228.611        
#> 5    12.4 3228.613        
#> 6    16.5      Inf        
#> 7    20.6 3228.570        
#> 8    24.7 3228.570        
#> 9    28.8      Inf        
#> 10   32.9      Inf        
#> 11   37.0      Inf        
summary(f) 
#>  Genealogical placement GMRF model fit 
#> 
#> Phylogenetic tree with 100 tips and 98 internal nodes.
#> 
#> Tip labels:
#>   t18, t1, t16, t9, t7, t11, ...
#> 
#> Rooted; includes branch length(s).
#> Range of coefficients: 
#> -0.734776592154461 1.29018258506584
#> Precision parameter (log tau): 0.0999999999999996 
#> 
#>   logprecision    RMSCLO   Neff
#> 1          0.1 0.6358537 91.454
#> 
if (FALSE) { # \dontrun{
plot(f)
} # }