# Scott Merrill
# Treatment means
# Dec 3

rm(list = ls(all = TRUE))


##############################################################
# categorical data
####


# number.of.total.samples = 100
# Matricies have a problem with names and numbers. 
#    Data frames do not seem to. 
large.woody.debris = NULL
harvest.type=NULL
sim.data = data.frame (large.woody.debris = large.woody.debris, 
	harvest.type=harvest.type)

sim.data

# Harvest type suggests different distributions of 
#	remaining large woody debris

harvest.names = c("slash","burn","hack",
		"beetle colony","green logging")
harvest.means = c(2,3,3,5,7)
harvest.sd = c(1,1,1,1.5,2)
harvest.samples = c(15,25,20,23,17)

# the error above is variance within each harvest type.
# also the possibility of additional error such as measurement error

count = 1
for (x in 1:length(harvest.means)) { 
	for (y in 1:harvest.samples[x]) {
		sim.data[count,1] = (rnorm(mean=harvest.means[x],
			sd=harvest.sd[x],n=1)
			# still want to add in measurement error
			+ rnorm(mean = 0, sd = 1, n = 1))
		sim.data[count,2] = harvest.names[x]		
		count = count + 1
	}
}
hist(sim.data[,1])
sim.data


mod1 = aov(sim.data[,1] ~ as.factor(sim.data[,2]))
mod1
summary(mod1)
coefficients(mod1)

summary.lm(mod1) # looks at first alphabetical against others (beetle kill vs others)

# Suggested for looking at differences between treatments
TukeyHSD(mod1)





contrasts(as.factor(sim.data[,2]))

#########
## contrasts are in aov also
#### 


## # compares "burn" to all others?
contrast3=c(0,1,0,0,0)
contrast3 = cbind(contrast3,c(0,1,0,0,0))
contrast3 = cbind(contrast3,c(0,1,0,0,0))
contrast3 = cbind(contrast3,c(0,1,0,0,0))

contrasts(HT)=contrast3

mod4 = aov(sim.data[,1] ~ HT,contrasts=contrast3)
coefficients(mod4)
summary.lm(mod4)
contrasts(HT)

4.2-1.3

mean(sim.data[c(1:15,41:100),1])
mean(sim.data[c(16:40),1])

# TukeyHSD {stats}

require(graphics)


summary(mod1) # fm1 <- aov(breaks ~ wool + tension, data = warpbreaks))
TukeyHSD(mod1,"HT", ordered = TRUE)  #     "tension", ordered = TRUE)
plot(TukeyHSD(fm1, "tension"))


example(TukeyHSD)