require(diversitree) mytree = read.nexus("BEAST_Rosales_final.nex") nopoly = multi2di(mytree, random=TRUE) chardata = read.csv("Musse_HSHR.csv", header=TRUE) genusnames = chardata$Genus states = data.frame(chardata$D, chardata$W, chardata$R, chardata$S, row.names=genusnames) names(states) = c("D", "W", "R", "S") pruned = drop.tip(nopoly, "Pseudolmedea") pruned = drop.tip(pruned, "Zizyphus") lik = make.musse.multitrait(pruned, states) startp = starting.point.musse.multitrait(pruned, lik) lik.l = constrain(lik, lambdaD ~ lambda0, lambdaW ~ lambda0, lambdaR ~ lambda0, lambdaS ~ lambda0, lambdaDW ~ lambda0, lambdaDR ~ lambda0, lambdaDS ~ lambda0, lambdaWR ~ lambda0, lambdaWS ~ lambda0, lambdaRS ~ lambda0, lambdaDWR ~ lambda0, lambdaDWS ~ lambda0, lambdaDRS ~ lambda0, lambdaWRS ~ lambda0, lambdaDWRS ~ lambda0) lik.m = constrain(lik, muD ~ mu0, muW ~ mu0, muR ~ mu0, muS ~ mu0, muDW ~ mu0, muDR ~ mu0, muDS ~ mu0, muWR ~ mu0, muWS ~ mu0, muRS ~ mu0, muDWR ~ mu0, muDWS ~ mu0, muDRS ~ mu0, muWRS ~ mu0, muDWRS ~ mu0) lik.qD = constrain(lik, qD01.0 ~ qD01.W, qD01.0 ~ qD01.R, qD01.0 ~ qD01.S, qD01.0 ~ qD01.WR, qD01.0 ~ qD01.WS, qD01.0 ~ qD01.RS, qD01.0 ~ qD01.WRS, qD01.0 ~ qD10.W, qD01.0 ~ qD10.R, qD01.0 ~ qD10.S, qD01.0 ~ qD10.WR, qD01.0 ~ qD10.WS, qD01.0 ~ qD10.RS, qD01.0 ~ qD10.WRS) lik.qW = constrain(lik, qW01.0 ~ qW01.D, qW01.0 ~ qW01.R, qW01.0 ~ qW01.S, qW01.0 ~ qW01.DR, qW01.0 ~ qW01.DS, qW01.0 ~ qW01.RS, qW01.0 ~ qW01.DRS, qW01.0 ~ qW10.D, qW01.0 ~ qW10.R, qW01.0 ~ qW10.S, qW01.0 ~ qW10.DR, qW01.0 ~ qW10.DS, qW01.0 ~ qW10.RS, qW01.0 ~ qW10.DRS) lik.qR = constrain(lik, qR01.0 ~ qR01.D, qR01.0 ~ qR01.W, qR01.0 ~ qR01.S, qR01.0 ~ qR01.DW, qR01.0 ~ qR01.DS, qR01.0 ~ qR01.WS, qR01.0 ~ qR01.DWS, qR01.0 ~ qR10.D, qR01.0 ~ qR10.W, qR01.0 ~ qR10.S, qR01.0 ~ qR10.DW, qR01.0 ~ qR10.DS, qR01.0 ~ qR10.WS, qR01.0 ~ qR10.DWS) lik.qS = constrain(lik, qS01.0 ~ qS01.D, qS01.0 ~ qS01.W, qS01.0 ~ qS01.R, qS01.0 ~ qS01.DW, qS01.0 ~ qS01.DR, qS01.0 ~ qS01.WR, qS01.0 ~ qS01.DWR, qS01.0 ~ qS10.D, qS01.0 ~ qS10.W, qS01.0 ~ qS10.R, qS01.0 ~ qS10.DW, qS01.0 ~ qS10.DR, qS01.0 ~ qS10.WR, qS01.0 ~ qS10.DWR) lik.qall = constrain(lik, qD01.0 ~ qD01.W, qD01.0 ~ qD01.R, qD01.0 ~ qD01.S, qD01.0 ~ qD01.WR, qD01.0 ~ qD01.WS, qD01.0 ~ qD01.RS, qD01.0 ~ qD01.WRS, qD01.0 ~ qD10.W, qD01.0 ~ qD10.R, qD01.0 ~ qD10.S, qD01.0 ~ qD10.WR, qD01.0 ~ qD10.WS, qD01.0 ~ qD10.RS, qD01.0 ~ qD10.WRS, qR01.0 ~ qR01.D, qR01.0 ~ qR01.W, qR01.0 ~ qR01.S, qR01.0 ~ qR01.DW, qR01.0 ~ qR01.DS, qR01.0 ~ qR01.WS, qR01.0 ~ qR01.DWS, qR01.0 ~ qR10.D, qR01.0 ~ qR10.W, qR01.0 ~ qR10.S, qR01.0 ~ qR10.DW, qR01.0 ~ qR10.DS, qR01.0 ~ qR10.WS, qR01.0 ~ qR10.DWS, qR01.0 ~ qR01.D, qR01.0 ~ qR01.W, qR01.0 ~ qR01.S, qR01.0 ~ qR01.DW, qR01.0 ~ qR01.DS, qR01.0 ~ qR01.WS, qR01.0 ~ qR01.DWS, qR01.0 ~ qR10.D, qR01.0 ~ qR10.W, qR01.0 ~ qR10.S, qR01.0 ~ qR10.DW, qR01.0 ~ qR10.DS, qR01.0 ~ qR10.WS, qR01.0 ~ qR10.DWS, qS01.0 ~ qS01.D, qS01.0 ~ qS01.W, qS01.0 ~ qS01.R, qS01.0 ~ qS01.DW, qS01.0 ~ qS01.DR, qS01.0 ~ qS01.WR, qS01.0 ~ qS01.DWR, qS01.0 ~ qS10.D, qS01.0 ~ qS10.W, qS01.0 ~ qS10.R, qS01.0 ~ qS10.DW, qS01.0 ~ qS10.DR, qS01.0 ~ qS10.WR, qS01.0 ~ qS10.DWR) fit = find.mle(lik, startp[argnames(lik)]) fit.l = find.mle(lik.l, startp[argnames(lik.l)]) fit.m = find.mle(lik.m, startp[argnames(lik.m)]) fit.qD = find.mle(lik.qD, startp[argnames(lik.qD)]) fit.qW = find.mle(lik.qW, startp[argnames(lik.qW)]) fit.qR = find.mle(lik.qR, startp[argnames(lik.qR)]) fit.qS = find.mle(lik.qS, startp[argnames(lik.qS)]) fit.qall = find.mle(lik.qall, startp[argnames(lik.qall)]) fit summary(fit) anova(fit, equal.l=fit.l) anova(fit, equal.m=fit.m) anova(fit, equal.q=fit.qD) anova(fit, equal.q=fit.qW) anova(fit, equal.q=fit.qR) anova(fit, equal.q=fit.qS) anova(fit, equal.q=fit.l)