library(emdbook) # for the datasets
library(plotrix)
library(gplots) 


data(ReedfrogPred)
summary(ReedfrogPred)
head(ReedfrogPred,10)
ReedfrogPred[1:10,]
bwplot(propsurv~density|pred*size,data=ReedfrogPred,horizontal=F)
boxplot(propsurv~size*density*pred,data=ReedfrogPred)


#### Sorting Data
attach(ReedfrogPred)

propsurv

sort(propsurv,decreasing=T)

ReedfrogPred.sort=ReedfrogPred[order(propsurv,decreasing=T),]
head(ReedfrogPred.sort)
ReedfrogPred.sort[1:6,]
detach()


### Visualizing tricky data
data(SeedPred)
windows();par(mfrow=c(3,1))
plot(SeedPred$available,SeedPred$taken)
plot(taken~available,data=SeedPred)#lattice, roughly identical to above plot
plot(jitter(SeedPred$available),jitter(SeedPred$taken)) # Adds noise, can see relative density


### Contingency tables, data manipulation, and visualization
t1=table(SeedPred$available,SeedPred$taken)
t1
plot(t1)
balloonplot(t1)#gplots package
dotchart(t1)#gplots package
attach(SeedPred)
n.by.avail=table(available)
frac.taken=SeedPred$taken/SeedPred$available
hist(frac.taken)
mean.frac.by.avail=tapply(frac.taken,available,mean,na.rm=TRUE)



### Ashton's Stuff

setwd("C:/Users/Matthew/Downloads")
caribmang = read.csv(file="caribmang.csv", head=TRUE, sep=",") #read the full dataset
summary(caribmang)                                             #show data in window to verify
names(caribmang)
names.fish=names(caribmang)
j=1
individ=sample(7:39,6,replace=F)
windows()
op <- par(mfrow=c(2,2))                #specify 2 plots side by side - 2 rows keeps them square
for (i in individ){


fish.sub =caribmang[,c(2,6,i)]   #subset out the 2 env variables plus one fish
fish.lab = names.fish[i]             #save the fish name for labeling files
mod2=paste(fish.lab,"~MgPerCoast",sep="")
##    Create pdf file with 2 figures side-by-side
##    (1)Fish in L vs. S islands
##    (2)Linear fit of Fish by Mangrove Area per km of Coastline

ifelse(j==2,{windows();op <- par(mfrow=c(2,2));j<-1},j<-2)#specify 2 plots side by side - 2 rows keeps them square,filling with 2 plots
plot(fish.sub[,1],fish.sub[,3],xlab="Island Group", ylab=fish.lab)
linear = lm(mod2, data=fish.sub)
plot(fish.sub[,2],fish.sub[,3],xlab=names.fish[6], ylab=fish.lab)
abline(linear)
                               #reset the plot options
                            #turn off the device window so the pdf will finish

ttest = t.test(fish.sub[,3] ~ IsGroup, data=fish.sub)                            
cat("These are results for the fish: ", fish.lab, ", run on ", date(),fill=T)
cat("T=",ttest$statistic," on ",ttest$parameter," degrees of freedom, P=",ttest$p.value,fill=T,sep="")
cat("-----------------------------------------------------------------------------------------------",fill=T)
}


### Alternative:
### Suppose we just to know which ones where T was significant (although I don't want to encourage this!!!!)
output.data=c(NA,NA)
for(i in 7:39){
fish.sub =caribmang[,c(2,6,i)]   #subset out the 2 env variables plus one fish
fish.lab = names.fish[i]             #save the fish name for labeling files
parm= t.test(fish.sub[,3] ~ IsGroup, data=fish.sub)$p.value
if(parm<0.05)output.data=rbind(output.data,c(fish.lab,round(parm,5)))
}
output.data=output.data[-1,]
output.data
output.data[,2]=as.numeric(output.data[,2])



########## Logistic Regression
data(ReedfrogFuncresp)
attach(ReedfrogFuncresp)
Frac.killed=Killed/Initial
plot(Initial,Frac.killed)
min(Initial)
max(Initial)
seq.mod=seq(5,100,.5)
mod=exp(mod$coef[1]+mod$coef[2]*seq.mod)/(1+exp(mod$coef[1]+mod$coef[2]*seq.mod))
mod=glm(cbind(Killed,Initial-Killed)~Initial,family="binomial")
mod.plot=exp(mod$coef[1]+mod$coef[2]*seq.mod)/(1+exp(mod$coef[1]+mod$coef[2]*seq.mod))
lines(seq.mod,mod.plot)