###Exercicio 6


#1. - 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

B. Idem para uma distribuição normal com média 7.5 e desvio padrão 3.2


> dados = data.frame(c(rep(NA,10)),c(rep(NA,10)))
> 
> for(i in 1:1)
+ 	{
+ 	dados[,i] = rnorm(10,mean = 6, sd = 3)
+ 	dados[,i+1] = rnorm(10,mean = 7.5, sd = 3.2)
+ 	}
> names(dados) = c("A","B")
> dados
           A          B
1   1.651109 10.9499050
2   7.957656  0.2897166
3   7.436934  7.0553977
4   8.557766  7.2341411
5   2.518030  7.3013053
6   3.658450 13.5529526
7   4.184556  2.5290620
8   9.220436 12.6169861
9  -2.238932  3.6718468
10 13.705996  6.3459437



DICA: utilize a função rnorm()

#1.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!
   
> ##Medias das amostras diferentes são diferentes ?
> 
> source("simula.r")
> 
> hist(dados$A, freq=FALSE)
> hist(dados$B, freq=FALSE)
> mean.A = mean(dados$A)
> mean.B = mean(dados$B)
> 
> dif=abs(mean(dados$A)-mean(dados$B))
> dif
[1] 1.489526
> 
> dados.AB= as.vector(c(dados$A,dados$B))
> 
> hist(dados.AB, freq=FALSE)
> curve(exp=dnorm(x, mean=mean(dados.AB),sd=sd(dados.AB)), col="red", add=T)
Warning message:
In if (add) lines(x, y, type = type, ...) else plot(x, y, type = type,  :
  a condição tem comprimento > 1 e somente o primeiro elemento será usado
> 
> mean = mean(dados.AB)
> sd = sd(dados.AB)
> 
> dif.AB = simula(dados$A,dados$B,nsim=2000)
simula uma distribuição nula, no caso de teste unicaudal o primeiro vetor deve ser o dos dados que seriam maiores.
 As opçoes de teste são (entre aspas):
         uni (normal unicaudal)
         bi (normal bicaudal) 
         t (distribuição t)

 Diferença absoluta observada entre as médias das variáveis =  1.5 
Houve 50 ou mais avisos (use warnings() para ver os primeiros 50)
> table(dif.AB)
dif.AB
  0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 
 47  80  85  86  67  88  83  77  70  83  61  69  63  62  52  84  58  53  45  53 
  2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9   3 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 
 50  55  28  42  38  47  37  36  35  27  30  15  18  15  20  17   9  17  12   9 
  4 4.1 4.2 4.3 4.4 4.5 4.7 4.8 4.9   5 5.2 5.3 5.4 5.5 5.6 5.8 5.9   6 6.1 6.3 
  6   5   7  10   3   1   7   4   8   5   1   1   5   1   3   1   2   2   1   1 
6.4 6.6 6.9 
  1   1   1 
> 
> n.maior=sum(dif.AB>=dif)
> n.maior
[1] 927
> 
> n.maior/length(dif.AB)
[1] 0.4635


#1.2 - Teste agora que a média da segunda amostra é maior que a primeira.

> dif.AB.uni = simula(dados$A,dados$B, teste="uni", nsim=2000)
simula uma distribuição nula, no caso de teste unicaudal o primeiro vetor deve ser o dos dados que seriam maiores.
 As opçoes de teste são (entre aspas):
         uni (normal unicaudal)
         bi (normal bicaudal) 
         t (distribuição t)

 Diferença absoluta observada entre as médias das variáveis =  1.5 
Houve 50 ou mais avisos (use warnings() para ver os primeiros 50)
> table(dif.AB.uni)
dif.AB.uni
-7.4 -6.2 -5.7 -5.6 -5.5   -5 -4.9 -4.8 -4.7 -4.6 -4.5 -4.4 -4.3 -4.2 -4.1   -4 
   1    1    1    1    2    2    1    2    5    2    2    3    2    5    1    4 
-3.9 -3.8 -3.7 -3.6 -3.5 -3.4 -3.3 -3.2 -3.1   -3 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 
   7    5    5    5   14   23   13   10   13   14   15   12   15   10   20   19 
-2.3 -2.2 -2.1   -2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1   -1 -0.9 -0.8 
  16   25   28   22   29   28   33   29   39   25   30   34   31   27   43   31 
-0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1    0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8 
  38   29   50   33   38   38   34   36   45   44   45   41   40   38   45   39 
 0.9    1  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9    2  2.1  2.2  2.3  2.4 
  47   27   35   39   24   36   36   36   26   29   20   24   26   21   19   18 
 2.5  2.6  2.7  2.8  2.9    3  3.1  3.2  3.3  3.4  3.5  3.6  3.7  3.8  3.9    4 
  12   16   10   17   16   14    9    7    9    7   11    9   11    3    5    6 
 4.1  4.2  4.3  4.4  4.5  4.6  4.9    5  5.1  5.2  5.4  5.5  5.6  5.8  6.1  6.6 
   4    4    3    5    2    3    3    3    1    1    2    1    1    1    1    1 
 8.6 
   1 
> 
> n.maior=sum(round(dif.AB.uni,1)>=round(dif,1))
> n.menor=sum(round(dif.AB.uni,1)<=round(dif*-1,1))
> n.maior
[1] 454
> n.menor
[1] 484
> 
> (n.maior+n.menor)/length(dif.AB)
[1] 0.469



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

> t.test(dados$A,dados$B) 

        Welch Two Sample t-test

data:  dados$A and dados$B 
t = -0.7497, df = 17.93, p-value = 0.4632
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -5.664977  2.685926 
sample estimates:
mean of x mean of y 
 5.665200  7.154726 

#Teste t
t.test(dados$A,dados$B) 

#Teste t na unha

vA=var(dados$A)
vA
vB=var(dados$B)
vB
nA=length(dados$A)
nA
nB=length(dados$B)
nB

sAB=sqrt((vA/nA)+(vB/nB))
sAB
dif
tvalor=dif/sAB
tvalor

res.t= simula(dados$A,dados$B,nsim=2000, teste="t")

maior.menor.t=sum(res.t>=tvalor | res.t<=-tvalor)
maior.menor.t
prob.t=maior.menor.t/length(res.t)
prob.t

## compare com o resultado do teste t da função t.test()

t.test(dados$A,dados$B)

   
#1.4 - Não esqueça de fazer um gráfico para mostrar os dados…

plot(1:20,dados.AB, pch=rep(c(15,16),each=10),col=rep(1:2,each=10))
	for(i in 1:10)
	{
	lines(c(i,i),c(dados.AB[i],mean.A),col=1)
	}
	for(j in 11:20)
	{
	lines(c(j,j),c(dados.AB[j],mean.B),col=2)
	}

lines(c(1,10),c(mean.A,mean.A),col=1)
lines(c(11,20),c(mean.B,mean.B),col=2)



#2. - Utilizando os dados da planilha caixeta.csv

caixeta = read.csv("caixeta.csv")


#2.1 - Calcule os valores de área basal por fuste

str(caixeta)
head(caixeta)
caixeta$area.b = round(((((caixeta$cap/pi)^2)*pi)/4)*0.01,2)

   
#2.2 - Calcule os valores de área basal por amostra em cada uma das localidades

str(caixeta$parcela)
levels(caixeta$local)

caixeta.chauas = caixeta[caixeta$local=="chauas",]
caixeta.chauas.AB = tapply(caixeta.chauas$area.b, caixeta.chauas$parcela, sum)
caixeta.jureia = caixeta[caixeta$local=="jureia",]
caixeta.jureia.AB = tapply(caixeta.jureia$area.b, caixeta.jureia$parcela, sum)
caixeta.retiro = caixeta[caixeta$local=="retiro",]
caixeta.retiro.AB = tapply(caixeta.retiro$area.b, caixeta.retiro$parcela, sum)

caixeta.area.b = data.frame(caixeta.chauas.AB,caixeta.jureia.AB,caixeta.retiro.AB)
names(caixeta.area.b) = c("chauas","jureia","retiro")
caixeta.v = as.vector(c(caixeta.area.b$chauas,caixeta.area.b$jureia,caixeta.area.b$retiro))


mean.parcela = as.vector(apply(caixeta.area.b, 2, mean))
mean.var = as.vector(apply(caixeta.area.b, 2, var))

   

#2.3 - Produza gráficos para mostrar os dados

plot(1:15,caixeta.v, pch=rep(c(15,16,17),each=5),col=rep(1:3,each=5), ylab= "Area Basal (cm^2)", xlim=c(0,20))
	for(i in 1:5)
	{
	lines(c(i,i),c(caixeta.v[i],mean.parcela[1]),col=1)
	}
	for(j in 6:10)
	{
	lines(c(j,j),c(caixeta.v[j],mean.parcela[2]),col=2)
	}
	for(z in 11:15)
	{
	lines(c(z,z),c(caixeta.v[z],mean.parcela[3]),col=3)
	}

lines(c(1,5),c(mean.parcela[1],mean.parcela[1]),col=1)
lines(c(6,10),c(mean.parcela[2],mean.parcela[2]),col=2)
lines(c(11,15),c(mean.parcela[3],mean.parcela[3]),col=3)


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


caixeta.names = as.factor(c(rep("chauas",5),rep("jureia",5),rep("retiro",5)))
caixeta.v
vetor.cor = as.vector(rep(c("red","black","blue"), each=5))
vetor.cor[1]

#Grafico SS Total

plot(caixeta.v, pch=(rep(c(15,16,17),each=5)),col=vetor.cor, ylab="Area Basal", xlab="Observações")
	for(i in 1:15)
	{
	lines(c(i,i),c(caixeta.v[i],mean(caixeta.v)),col=vetor.cor[i])
	}
	abline(h=mean(caixeta.v))

## Calculo SS TOTAL

media.geral=mean(caixeta.v)
media.geral
dif.geral=caixeta.v-media.geral
dif.geral
sum(dif.geral)
round(sum(dif.geral),10)
ss.caixeta=dif.geral^2
ss.caixeta
ss.total=sum(ss.caixeta)
ss.total

##Somatorio Quadratico Intragrupos

vetor.media=c(rep(mean(caixeta.v[1:5]), each=5),rep(mean(caixeta.v[6:10]), each=5),rep(mean(caixeta.v[11:15]), each=5))

plot(caixeta.v, pch=(rep(c(15,16,17),each=5)),col=vetor.cor, main="Variação Intra Grupos",ylab="Variável Resposta", xlab="Observações")
	for(i in 1:15)
	{
	lines(c(i,i),c(vetor.media[i],caixeta.v[i]),col=vetor.cor[i])
	}
	lines(c(1,5),c(vetor.media[1],vetor.media[1]),col=vetor.cor[1])
	lines(c(6,10),c(vetor.media[6],vetor.media[6]),col=vetor.cor[6])
	lines(c(11,15),c(vetor.media[11],vetor.media[11]),col=vetor.cor[11])

caixeta.v
vetor.media
	for (j in 1:5)
	{
	ss.caixeta[j]=(caixeta.v[j]-vetor.media[j])^2
	ss.caixeta[j+5]=(caixeta.v[j+5]-vetor.media[j+5])^2
	ss.caixeta[j+10]=(caixeta.v[j+10]-vetor.media[j+10])^2
	}
	ss.caixeta
ss.intra = sum(ss.caixeta)
ss.intra


#Somatorio da Variacao Quadratica entre grupos

plot(vetor.media,pch=(rep(c(15,16,17),each=5)),col=vetor.cor, main="Variação Entre Grupos",ylab="Variável Resposta", xlab="Observações")
		for(i in 1:15)
		{
		lines(c(i,i),c(vetor.media[i],mean(caixeta.v)),col=vetor.cor[i])
		}
	abline(h=mean(vetor.media))
	
#### Cálculo dos valores

media.solos=apply(solos,2,mean)
media.solos
media.geral
ss.entre=10*sum((media.solos-media.geral)^2)

	
	ss.entre=c(rep(NA,15))
	for (i in 1:15)
	{
	ss.entre[i]=(vetor.media[i]-mean(caixeta.v))^2
	}
	ss.entre = sum(ss.entre)

ss.total = ss.intra+ss.entre
ss.total


## Calculo de F

gl.total = lenght(caixeta.v)-1
gl.entre = 3-1
gl.intra = gl.total-gl.entre


ms.entre=ss.entre/gl.entre
ms.intra=ss.intra/gl.intra
ms.entre
ms.intra

F.caixeta.dap=ms.entre/ms.intra
F.caixeta.dap

curve(expr=df(x, gl.entre,gl.intra),main="Disribuição F de Fisher", xlab="Valor F",ylab="Densidade",xlim=c(0,10))
abline(v=F.caixeta.dap,col="red")

##Calculo do valor de P

p.caixeta.dap=pf(F.caixeta.dap,gl.entre,gl.intra, lower.tail=FALSE)
p.caixeta.dap

#2.5 - Calcule a tabela de anova com a função aov().

caixeta.names = factor(c(rep("chauas",5),rep("jureia",5),rep("retiro",5)))

aov(caixeta.v~caixeta.names)

   
#2.6 - Qual é a porcentagem de variação explicada pela localidade nesse caso?

> p.explicada.local = ss.entre/ss.total*100
> p.explicada.local
[1] 34.39258