### Exercício 6 ####

####Crie seus dados#####
#¦Crie dois conjuntos de dados: 
#A. Dez observações de uma amostra de uma distribuição normal com média 6 e desvio padrão 3 

ca=rnorm(10,mean=6,sd=3)
ca
class(ca)

#B. Idem para uma distribuição normal com média 7.5 e desvio padrão 3.2 
cb=rnorm(10,mean=7.5, sd=3.2)
cb
class(cb)

plot(ca,cb, col=c("red","black"), pch=c(17,5), cex=1.5, ylab="freq", xlab="")
boxplot(ca,cb,ylab="freq", col="gray")
#DICA: utilize a função rnorm() 

#####fazendo igual ao tutorial######

obj.dados=c(ca,cb) 
obj.dados
fator=factor(rep(c("a","b"),each=10))
boxplot(obj.dados~fator)

#diferença entre médias das coletas de ca e cb
(mean(rnorm(10,mean=mean(ca),sd=sd(ca)))-(mean(rnorm(10,mean=mean(cb),sd=sd(cb)))))

resulta=rep(NA,10)
for (i in 1:1000)
	{
	resulta[i]= (mean(rnorm(10,mean=mean(ca),sd=sd(ca)))-(mean(rnorm(10,mean=mean(cb),sd=sd(cb)))))}

obs=mean(ca)-mean(cb)
obs 

p= sum(resulta>=obs)/length(resulta)
p# 0.494 - não posso afirmar que são diferentes, fico com a hip nula


## fazendo com o simula.r

#1.Utilize a função simula.r e teste a hipótese que as médias das amostras são diferentes. Não esqueça de usar a função source() para carregar a função!

source("http://ecologia.ib.usp.br/bie5782/lib/exe/fetch.php?media=bie5782:01_curso2009:simula.r")
simula# function(dados1,dados2, nsim=1000, teste="bi")

sim1=simula(ca,cb, nsim=1000, teste="bi")
summary(sim1)

p2=((sum(sim1>=obs))/length(sim1))
p2 #0.21 - fico com hipótese nula, minha chance de errar ao afirmar que ca é diferente de cb é muito alta (21%)

#2.Teste agora que a média da segunda amostra é maior que a primeira. Compare os valores obtidos. Atenção: nunca faça isso na vida real! Sua hipótese deve estar definida a priori, antes de iniciar o teste! 

simula2=simula(cb,ca,nsim=1000, teste="uni")
simula2

n.maior=sum((simula2)>=(obs))
n.maior
n.menor=sum((simula2)<=(obs*-1))
n.menor

p.uni=(n.menor+n.maior)/length(simula2)
p.uni # 0.229 - também não aceito h0

#3.Utilize agora a função t.test() para testar as mesmas hipóteses. Os resultados são iguais?

t.test(ca,cb)# p=0.2117 
#sim, continuo com a h0 !

#4.Você fez uma Análise Exploratória dos Dados antes de fazer o teste? DEVERIA! Acostume-se a criar intimidade com os dados antes de fazer qualquer teste. Como diria o professor Rodrigo Pereira: “Leve os dados para passear, tomem um cafezinho, tornem-se íntimos!”
#fiz o boxplot !

#5.Sempre faça um diagnóstico das premissas do teste! Quais são as premissas do teste na função simula? E na função t.test()? Faça o diagnóstico gráfico dessas premissas! 
# as duas tem as premissas de igualdade de variâncias e distribuição normal !

# normalidade

qqnorm(ca)
qqline(ca)
qqnorm(cb)
qqline(cb)

#variancia
bartlett.test(list(ca,cb))#p=0.2294 - são iguais

####Caixeta de NOVO?!######

#Utilizando os dados da planilha caixeta.csv: 

caix=file.choose()
caixeta=read.table(caix,header=TRUE, sep=",",as.is=TRUE)
head(caixeta)

# 1.Calcule o valor da área basal para cada indivíduo (somatória dos fustes). 

raio=(caixeta$cap)/(2*pi)
caixeta$bas=base::pi*(raio^2)


# 2.Calcule o valor da área basal por amostra em cada uma da localidades (somatória da área basal dos individuos por amostra) 
help(aggregate)


##aggregate###

basal.par=aggregate(caixeta$bas,list(caixeta$local,caixeta$parcela),sum)
# tem que colocar o local tb, caso contrário eu junto as parcelas dos locais diferentes !
# seria melhor fazer uma tabela com ns diferentes nas parcelas 
basal.par

###tapply###

basal.par2=tapply(caixeta$bas,list(caixeta$local,caixeta$parcela),sum)
basal.par2
# para ordenar do menor para o maior
help(sort)#sort(x, decreasing = FALSE, ...)
sort(basal.par[,3],decreasing = FALSE )# todas as linhas, coluna 3

#comparando
head(caixeta)
basal.par
class(basal.par)#data.frame
basal.par2
class(basal.par2)#matriz-prefiro esse formato para visualizar as parcelas !


# 3.Produza gráficos para mostrar os dados e os diferentes desvios relacionados às amostras e as localidades, sendo as localidades o seu fator de interesse (variação total, intra grupos e entre grupos)

par(mfrow=c(1,3))
basal.par2
boxplot(basal.par2[1,],ylmfrow=c(2,3),ylim=c(6000,1e+06), xlab="chauas")
boxplot(basal.par2[2,], xlab="jureia")
boxplot(basal.par2[3,],xlab="retiro")

#tentando usar matpot
cores=c("black","black","black","black","black","red","red","red","red","red","blue","blue","blue","blue","blue")
matplot(basal.par2, las=TRUE, col=cores)# nao deu mt certo =/


4.Calcule os valores de uma tabela Anova para esses dados 
#sendo a variável dependente a área basal/parcela e o tratamento 
#as localidades. Cada observação refere-se a uma amostra ou parcela.

#Desvios Quadráticos 

#separando dados dos locais 
basal.par
chauas=basal.par[c(1,4,7,10,13),3]
chauas
retiro=basal.par[c(3,6,9,12,15),3]
retiro
jureia=basal.par[c(2,5,8,11,14),3]
jureia

data=data.frame(chauas,retiro,jureia)
data

#medias de cada local/ de cada grupo
media.cha= mean(chauas)
media.ju=mean(jureia)
media.re=mean(retiro)

#ou

medias.loc=mean(data)
medias.loc

#media geral/ media dos grupos

media.geral=mean(medias.loc)
media.geral


# SS TOTAL
ss=medias.loc-media.geral
ss
sum(ss)
round(sum(ss),5)
SS.loc=ss^2
SS.loc
SS.total=sum(SS.loc)
SS.total

# SS INTRAGRUPOS
#pontos de cada grupo - media dentro do grupo
SS.chauas=sum((chauas-media.cha)^2)
SS.chauas
SS.retiro=sum((retiro-media.re)^2)
SS.retiro
SS.jureia=sum((jureia-media.ju)^2)
SS.jureia
SS.intra=SS.chauas+SS.retiro+SS.jureia
SS.intra

# SS ENTREGRUPOS
#media do grupo- media geral, multiplicado pelo tamanho do grupo e elevado ao quadrado
# todos tem 5 valores

grupos= (5*(medias.loc-media.geral)^2)
SS.entre=sum(grupos)
SS.entre

#  Estatistica F 
MS.entre=SS.entre/2
MS.intra=SS.intra/12
MS.entre
MS.intra
F=MS.entre/MS.intra
F
pf(F, 2, 12, lower.tail=FALSE) #0.07974

# AOV

help(aov)#aov(formula, data = NULL, projections = FALSE, qr = TRUE, contrasts = NULL, ...)
head(basal.par)

locais=factor(basal.par$Group.1)
class(locais)
locais
a.basais=basal.par$x
a.basais

novop= aov(a.basais~locais)
novop
summary(novop) # P=0.07975 ! =D foi igual !!!!


6.Qual é a porcentagem de variação explicada pela localidade nesse caso?

# é proporcional ao p =0.079 - explica quase 80%