source("mleFunctions.r")
source("bbmle.R")
data <- read.csv("sokalAndRolhfRatData.csv")

ll <- function (m,s)
-sum(dnorm(data$consumed,mean=m,sd=s,log = TRUE))
guess <- list(m = 6, s = 1)
fit <- mle(ll, start=guess, method="Nelder-Mead")

summary(fit)
mean(data$consumed)
#note: mean is off, try random restart
fit <- mle.random(ll, guess, lower.guess = c(min(data$consumed),0), upper.guess=c(max(data$consumed),100))
summary(fit)
comp <- compare(fit, name="null")

#categorical independent variable
#be sure vector in same order as this
levels(data$status)
ll <- function (m.f,m.r,s) with(data,{
    m = c(m.f,m.r);
    -sum(dnorm(consumed,mean=m[status],sd=s,log = TRUE))
    })
guess <- list(m.f = 500, m.r = 450, s = 80)
fit <- mle(ll, start=guess, method="Nelder-Mead")
comp <- compare(fit, name="status main", prev=comp)
#
levels(data$sex)
ll <- function (m.female,m.male,s) with(data,{
    m = c(m.female,m.male);
    -sum(dnorm(consumed,mean=m[sex],sd=s,log = TRUE))
    })
guess <- list(m.female = 500, m.male= 450, s = 80)
fit <- mle(ll, start=guess, method="Nelder-Mead")
comp <- compare(fit, name="sex main", prev=comp)
#
with(data,levels(sex:status))
ll <- function (f.f,f.r,m.f,m.r,s) with(data,{
    m = c(f.f,f.r,m.f,m.r);
    -sum(dnorm(consumed,mean=m[sex:status],sd=s,log = TRUE))
    })
guess <- list(f.f = 500, f.r= 450, m.f=500, m.r=500, s = 80)
fit <- mle(ll, start=guess, method="Nelder-Mead")
comp <- compare(fit, name="sex:main", prev=comp)

compare.likRatio(comp,1,2)
compare.likRatio(comp,2,3)

## slight automation, at least you don't have to enter each starting condition individually
with(data,levels(sex:status))
ll <- function (f.f,f.r,m.f,m.r,s) with(data,{
    m = c(f.f,f.r,m.f,m.r);
    -sum(dnorm(consumed,mean=m[sex:status],sd=s,log = TRUE))
    })
interact <- with(data,interaction(sex,status))
guess <- rep(500,nlevels(interact))
names(guess) <- levels(interact)
guess <- c(as.list(guess),s=80)
fit <- mle(ll, start=guess, method="Nelder-Mead")

## automation notes
#this can be used to set the list of arguments expected in ll
formals(ll) <- alist(a=)