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)