```
library(ape)
```

```
## Warning: package 'ape' was built under R version 2.15.1
```

```
lcx <- 1 ## contour lab cex
dumpfn <- "acedump.RData"
if (file.exists(dumpfn)) load(dumpfn)
mytree <- read.nexus("monkey_tree.nex")
copy_num <- c(35, 24, 6)
names(copy_num) <- c("human", "chimp", "rhesus")
source("ace.R") ## hacked versions of curve3d, apply2d, and ace
```

```
## Loading required package: plyr
```

```
max_copy_num = 80
```

Set up the transition matrix:

```
copy_num <- factor(copy_num, levels = 1:max_copy_num)
transitions <- matrix(0, max_copy_num, max_copy_num)
transitions[(row(transitions) - col(transitions)) == -1] <- 1
transitions[(row(transitions) - col(transitions)) == 1] <- 2
```

Regular `ace`

fails:

```
out <- ace(copy_num, mytree, type = "discrete", model = transitions)
```

```
## Error: non-finite value supplied by 'nlm'
```

Now we'll use the hacked version to return just the deviance function so we can play with it:

```
info <- ace2(copy_num, mytree, type = "discrete", model = transitions, devInfoOnly = TRUE)
```

See what's in there:

```
names(info)
```

```
## [1] "np" "ip" "fn"
```

```
info$np
```

```
## [1] 2
```

```
info$ip
```

```
## [1] 0.1
```

Generate deviance surface (by brute force):

```
if (!exists("cc")) {
cc <- curve3d(info$fn(c(x, y)), xlim = c(0.01, 100), ylim = c(0.01, 100),
log = "xy", sys3d = "none")
}
```

Picture:

```
with(cc, image(x, y, z, log = "xy"))
with(cc, contour(x, y, z, log = "xy", add = TRUE))
```

Things are so horrid that the poor thing can't build a contour plot effectively.

Just for comparison we'll try `persp`

— since it can't handle infinite values

we'll replace them:

```
maxz <- with(cc, ifelse(z == Inf, max(is.finite(z)), z))
with(cc, persp(x, y, maxz))
```

Still pretty ugly!

… but there are some reasonable values in the results:

```
hist(na.omit(c(cc$z)), col = "gray", main = "")
```

Try again, over a narrower range:

```
if (!exists("cc2")) {
cc2 <- curve3d(info$fn(c(x, y)), xlim = c(0.1, 10), ylim = c(0.1, 10), log = "xy",
sys3d = "none")
}
```

Better, but still ugly:

```
with(cc2, image(x, y, z, log = "xy"))
with(cc2, contour(x, y, z, log = "xy", labcex = lcx, add = TRUE))
```

If we try Nelder-Mead (a more robust algorithm), and start within

the OK region, we can actually get a reasonable answer:

```
(opt1 <- optim(fn = info$fn, par = c(1, 1), method = "Nelder-Mead"))
```

```
## $par
## [1] 68.83 74.97
##
## $value
## [1] 15.23
##
## $counts
## function gradient
## 83 NA
##
## $convergence
## [1] 0
##
## $message
## NULL
```

Now let's check the vicinity of those results:

```
if (!exists("cc3")) {
cc3 <- curve3d(info$fn(c(x, y)), xlim = c(50, 80), ylim = c(50, 80), sys3d = "none")
}
```

```
with(cc3, image(x, y, z))
with(cc3, contour(x, y, z, add = TRUE, labcex = lcx))
with(cc3, contour(x, y, z, add = TRUE, levels = 15:20, lty = 2, labcex = lcx))
points(opt1$par[1], opt1$par[2], pch = 16)
```

If all we care about is the MLE, then we're done. I thought

we might want to know something about the confidence intervals, though \ldots

We can calculate the local curvature (Hessian) and from

that get the correlation matrix of the parameters:

```
library(numDeriv)
h1 <- hessian(info$fn, opt1$par)
v <- solve(h1)
cov2cor(v)
```

```
## [,1] [,2]
## [1,] 1.000 0.998
## [2,] 0.998 1.000
```

\( x \) and \( y \) are almost perfectly correlated.

We would be better off reparameterizing this in terms of

(either) \( x+y \) and \( x-y \) (I also considered \( xy \) and \( x/y \),

but the former seems to work OK).

```
fn2 <- function(sumxy, diffxy) {
x <- (sumxy + diffxy)/2
y <- sumxy - x
info$fn(c(x, y))
}
```

Checking:

```
info$fn(c(1, 1))
```

```
## [1] 56.97
```

```
fn2(2, 0)
```

```
## [1] 56.97
```

```
info$fn(c(1, 0.7))
```

```
## [1] 61.52
```

```
fn2(1.7, 0.3)
```

```
## [1] 61.52
```

I drew a few surfaces, tweaking the limits outward each time:

```
if (!exists("cc4")) {
cc4 <- curve3d(fn2(x, y), xlim = c(20, 400), ylim = c(-20, 5), sys3d = "none")
}
```

In order to fit the function on this scale we need to make a little wrapper:

```
fn3 <- function(p) {
fn2(p[1], p[2])
}
opt3 <- optim(fn = fn3, par = c(150, -5), method = "Nelder-Mead")
```

A picture:

```
with(cc4, image(x, y, z, xlab = "x+y", ylab = "x-y", zlim = c(15, 25)))
with(cc4, contour(x, y, z, add = TRUE, labcex = lcx, levels = 14:25))
with(cc4, contour(x, y, z, add = TRUE, labcex = lcx, levels = seq(15, 16, by = 0.1),
lty = 2))
points(opt3$par[1], opt3$par[2], pch = 16)
```

Expand the range a bit more: to get 95% bivariate confidence intervals

we have to get up to (min+6) \( \approx \) 21.3 …

```
if (!exists("cc5")) {
cc5 <- curve3d(fn2(x, y), xlim = c(20, 600), ylim = c(-25, 5), sys3d = "none")
}
```

```
with(cc5, image(x, y, z, xlab = "x+y", ylab = "x-y", zlim = c(15, 25)))
with(cc5, contour(x, y, z, add = TRUE, labcex = lcx, levels = 14:25))
with(cc5, contour(x, y, z, add = TRUE, labcex = lcx, levels = seq(15, 16, by = 0.1),
lty = 2))
```

We might try reparameterizing this as \( \{1/(x+y),x-y\} \), so that we can try the boundary as \( x+y \to \infty \)?

```
fn2(1000, -10) ## still not above 16!
```

```
## [1] 15.92
```

Compute cross-sections of the deviance surface for a couple of large values of \( x+y \):

```
xydiffvec <- seq(-200, 5, by = 1)
if (!exists("bigslice")) bigslice <- sapply(xydiffvec, function(d) fn2(5000,
d))
plot(xydiffvec, bigslice, type = "l")
```

The midpoint of the surface (i.e. the best value of \( x-y \) for a given value

of \( x+y \)) keeps decreasing as we go out:

```
xydiffvec2 <- seq(-500, 0, by = 2)
if (!exists("bigslice2")) bigslice2 <- sapply(xydiffvec2, function(d) fn2(10000,
d))
plot(xydiffvec2, bigslice2, type = "l")
```

Save lengthy computation results:

```
if (!file.exists(dumpfn)) save(list = ls(pattern = "^(cc|big)"), file = dumpfn)
```