# ------------ #
# Lecture 8  
# ------------ #

bacteria.table <- read.table("Bacteria.txt", header=T)
attach(bacteria.table)

##PLOT: bacteria vs. time -- exponential decay?
#jpeg('bacteria.jpg', height=600, width=600)
plot(bacteria.table, pch=23, cex=2, bg='orange')
#dev.off()
##PLOTEND

# Fit model with untransformed data

bacteria.lm <- lm(N_t ~ t)
plot(bacteria.lm, cex=2, pch=23, bg='orange')
plot(bacteria.table, pch=23, cex=2, bg='orange')
lines(t, fitted(bacteria.lm), lwd=2, col='red')

# Fit model with log-transformed data

bacteria.log.lm <- lm(log(N_t) ~ t)
par(mfrow=c(2,2))
plot(bacteria.log.lm, cex=2, pch=23, bg='orange')


par(mfrow=c(1,1))
plot(bacteria.table, pch=23, cex=2, bg='orange')
lines(t, fitted(bacteria.lm), lwd=2, col='red')
lines(t, exp(fitted(bacteria.log.lm)), lwd=2, col='green')


# ----- Managers data ----    

manager = read.table("Manager.txt",header=T)
attach(manager)

manager.lm=lm(Y~X)

#jpeg("manager_fit.jpg",width=600,height=600)
plot(X,Y,col='blue')
abline(manager.lm)
#dev.off()

plot(X,rstandard(manager.lm))

plot(X,abs(resid(manager.lm)))

# Do Brown-Forsythe Test.


BrownForsythe <- function(r,Group1,Group2){
    n1 = length(Group1)
    n2 = length(Group2)
    r1 = r[Group1]
    r2 = r[Group2]
    rmed1 = median(r1)
    rmed2 = median(r2)
    D1 <- sum( abs( r1 - rmed1 )) / n1
    D2 <- sum( abs( r2 - rmed2 )) / n2
    s <- sqrt( ( sum( ( abs(r1 - rmed1) - D1 )^2 ) + sum( ( abs(r2 - rmed2) - D2 )^2 ) ) / (n1+n2-2) )
    t <- abs(( D1 - D2 ) / ( s * sqrt( 1/n1 + 1/n2 ) ))
    pval = pt(t,n1+n2-2,lower.tail=FALSE)
    list(t=t,pval=pval)
}

Group1 = 1:floor(length(X)/2)
Group2 = (floor(length(X)/2)+1):length(X)
BF = BrownForsythe(resid(manager.lm), Group1,Group2)

low=lowess(X, abs(resid(manager.lm)))
lines(low,col="blue")

X2=1/X;
Y2=Y/X;
manager.lm2 = lm(Y2~X2)

plot(X2,Y2)
abline(manager.lm2)

plot(X,abs(resid(manager.lm2)))
BF2 = BrownForsythe(resid(manager.lm2), Group1,Group2)

Y3 = log(Y)
manager.lm3 = lm(Y3~X)

plot(X,Y3)
abline(manager.lm3)
plot(X,abs(resid(manager.lm3)))

X3 = X^2
manager.lm4 = lm(Y3~X+X3)
plot(X,Y3)
points(X,fitted(manager.lm4),pch=15,col="red")
plot(X,abs(resid(manager.lm4)))

# ---- Brain Data ----

brain=read.table("Brain.txt",header=T,sep='\t')
attach(brain)

plot(BodyWeight,BrainWeight)

# why don't we simply treat the outlying points as outliers?
plot(BodyWeight,BrainWeight,ylim=c(0 ,1000),xlim=c(0 ,1000))

#look at the distribution of data, note that they are 
# highly asymmetrical.
hist(BrainWeight,100)
hist(BodyWeight,100)
hist(log(BodyWeight),100)
hist(log(BrainWeight),100)

plot(log(BodyWeight),log(BrainWeight))
brain.log.lm = lm(log(BrainWeight)~log(BodyWeight))
abline(brain.log.lm)
dino=c(6,16,25)
points(log(BodyWeight[dino]),log(BrainWeight[dino]),col="red",pch=15)

plot(log(BodyWeight),log(BrainWeight))
for( i in c(1:28)){ text(log(BodyWeight[i]),log(BrainWeight[i]),Name[i]) }

plot(brain.log.lm)

throw=c(16,6,25)
brain.log.lm2 = update(brain.log.lm,subset=-throw)

plot(log(BodyWeight),log(BrainWeight))
abline(brain.log.lm)
abline(brain.log.lm2,col='red')
points(log(BodyWeight[dino]),log(BrainWeight[dino]),col="red",pch=15)
plot(brain.log.lm2)

# Plutonium Measurements

plutonium = read.table("Plutonium.txt",header=T,sep="\t")
plot(PlutoniumActivity, AlphaCountRate,pch=15)
throw = c(24)
plutonium.lm=lm(AlphaCountRate~PlutoniumActivity)

r = residuals(plutonium.lm)
Group1 = which(PlutoniumActivity<=5)
Group2 = which(PlutoniumActivity>5)
BM = BrownForsythe(r,Group1,Group2)

plutonium = plutonium[1:23,]
attach(plutonium)
plutonium.lm=lm(AlphaCountRate~PlutoniumActivity)
r = residuals(plutonium.lm)
Group1 = which(PlutoniumActivity<=5)
Group2 = which(PlutoniumActivity>5)
BM = BrownForsythe(r,Group1,Group2)


plot(PlutoniumActivity,sqrt(AlphaCountRate))
plutonium.lm2 = lm(sqrt(AlphaCountRate)~PlutoniumActivity)
plot(plutonium.lm2,pch=15)


plot(sqrt(PlutoniumActivity),sqrt(AlphaCountRate))
plutonium.lm3 = lm(sqrt(AlphaCountRate)~PlutoniumActivity)
plot(plutonium.lm3,pch=15)