library(MASS)
data(Animals) # load the data set animals from MASS
attach(Animals)
pdf("brain.pdf")
plot(brain ~ body, data=Animals, xlab="Body Weight (kg)", ylab="Brain Weight (g)", main="Original Units")
dev.off()
identify(Animals$body, Animals$brain, rownames(Animals))



out = identify(body, brain, rownames(Animals))
# save plot via file menu



brains.lm = lm(brain ~ body, data=Animals)
pdf("brains-resid.pdf")
par(mfrow=c(2,2))
plot(brains.lm)
dev.off()

row.names(Animals) == "Brachiosaurus"  # case 26

plot(lm(brain ~ body, data=Animals, subset= c(-26)))
plot(lm(brain ~ body, data=Animals, subset= c(-26, -6)))
plot(lm(brain ~ body, data=Animals, subset= c(-26, -6, -16)))

# no more dinosaurs, but now elephants are  the problems
pdf("brains-BC.pdf")
boxcox(brains.lm)
dev.off()

pdf("brain-resid-logY.pdf")
par(mfrow=c(2,2))
plot(lm(log(brain) ~ body, data=Animals))
dev.off()

pdf("brains-tran.pdf")
par(mfrow=c(2,2))
logbrain.lm = lm(log(brain) ~ log(body), data=Animals)
plot(logbrain.lm)
dev.off()

#use R^2 to compare the two (as long as same response)

pt(-max(abs(rstudent(logbrain.lm)))), 26)/2

#compare to Bonferroni adjusted alpha
.05/28

#outlier test fails - masking!

pdf("brain-log.pdf")
plot(brain ~ body, data=Animals, xlab="Body Weight (kg)", ylab="Brain Weight (g)", log="xy", main="Logarithmic Scale")
text(body[out], brain[out], row.names(Animals)[out])
dev.off()


logbrains.nodino.lm = lm(log(brain) ~ log(body) + I(row.names(Animals) == "Triceratops") +  I(row.names(Animals) == "Brachiosaurus") + I(row.names(Animals) == "Dipliodocus"), data=Animals)

library(xtable)
xtable(
anova(logbrains.nodino.lm, logbrain.lm)
)
xtable(summary(logbrains.nodino.lm))


library(BAS)
brains.bas = bas.lm(log(brain) ~ log(body) + diag(28), data=Animals, prior="hyper-g-n", a=3, modelprior=beta.binomial(1,28), method="MCMC", n.models=2^17, MCMC.it=2^18)
# check for convergence
pdf("bas-animals-conv.pdf")
plot(brains.bas$probne0, brains.bas$probs.MCMC)
dev.off()
pdf("brains-image.pdf")
image(brains.bas)
dev.off()

rownames(Animals)[c(6, 14, 16, 26)]


######## Nonlinear Regression ##########

x = c(2,4,6,8,10,24,28, 32)
y = c(1.63, 1.01, .73, .55, .41, .01, .06, .02)


conc.lm = lm(I(log(y) - log(30)) ~ x)

vhat = exp(-coef(conc.lm)[1])
khat = -coef(conc.lm)[2]

pdf("nonlinear.pdf", width=11,height=8.5)
par(mfrow=c(1,2))
plot(x, y)
lines(x, (30/vhat)*exp(-khat*x))

plot(fitted(conc.lm),residuals(conc.lm))
abline(h=0)
dev.off()

df = data.frame(y=y, x=x)
logconc.nlm = nls( log(y) ~ log((30/V)*exp(-k*x)), data=df, start=list(V=vhat, k=khat))
summary(logconc.nlm)

conc.nlm = nls( y ~ (30/V)*exp(-k*x), data=df, start=list(V=vhat, k=khat))
summary(conc.nlm)

pdf("nonlinear.pdf", width=11,height=8.5)
par(mfrow=c(1,2))
plot(x, y)
lines(x, exp(predict(logconc.nlm)), col=2, lty=2, lwd=1.5)
lines(x, predict(conc.nlm), col=4, lty=4, lwd=1.5)
legend(15, 1.5, legend=c("Log-Normal", "Normal"), col=c(2,4), lty=c(2,4 ))


plot(exp(predict(logconc.nlm)), y - exp(predict(logconc.nlm)), col=2, xlim=range(0, 1.5), ylim=c(-.10, .30), pch=15)
points(predict(conc.nlm),y - predict(conc.nlm), col=4, pch=16)
abline(h=0)
legend(0.01, .25, legend=c("Log-Normal", "Normal"), col=c(2,4), pch=c(15,16 ))
dev.off()

library(R2jags)



nlmodel= function() {

  v ~ dt(0, phi, 1)%_% T(0,)
  k ~ dt(0,phi, 1)%_% T(0,)
  #    b ~ dwish(I, 1)
  #    R ~ dwish(I, 1)
  #    lambda ~ dgamma(1/2, 1/2)
  for (i in 1:n) {
    y[i] ~ dnorm(mu[i], phi)
    mu[i] <- D/v * exp(- k * x[i])
  }
  phi <- pow(sigma, -2)
  sigma ~ dt(0.00000E+00, 1, 1) %_% T(0.00000E+00, )
  cl <- v*k
  t.5 <- log(2)/k
}




data = list(y=y, x=x,  n=length(y), D=30)

out = jags(model=nlmodel, data=data, param=c("v","k","cl","t.5", "sigma"), n.iter=30000, n.burnin=5000)

pairs(out$BUGSoutput$sims.matrix)
out   # summary

acf(out$BUGSoutput$sims.matrix[,"v"])