LUDMILA MARIA RATTIS TEIXEIRA

Exercícios Aula 6 - Teste de Hipótese e Simulação

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

DICA: utilize a função rnorm() 

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!

R:
> data.a=abs(round(rnorm(10,mean=6,sd=3))) 
> data.a
 [1]  4  7  9  4  3 14  2  5  6  7
> data.b=abs(round(rnorm(10,mean=7.5,sd=3.2)))
> data.b
 [1] 10  9 14  6  6  5  4  4  3  1
> media.data.a = mean(data.a) 
> media.data.a
[1] 6.1
> media.data.b = mean(data.b) 
> media.data.b
[1] 6.2
> conj. = c(data.a,data.b) 
> conj. 
 [1]  4  7  9  4  3 14  2  5  6  7 10  9 14  6  6  5  4  4  3  1
> media.conj. = mean(conj.) 
> media.conj.
[1] 6.15
> sd.conj.= sd(conj.) 
> sd.conj.
[1] 3.558163
> dif.entre = abs(media.data.a-media.data.b) 
> dif.entre
[1] 0.1
> ##testando pelo boxplot:
> nomes.conj.=factor(rep(c("data.a","data.b"),each=10))
> nomes.conj.
 [1] data.a data.a data.a data.a data.a data.a data.a data.a data.a data.a
[11] data.b data.b data.b data.b data.b data.b data.b data.b data.b data.b
Levels: data.a data.b
> boxplot(conj.~ nomes.conj.)
> ##testando pelo simula.r( )
> resposta= rep(NA, 10)
> for(i in 1:10)
+ {
+ cat("\n\t", i) 
+ } 

         1
         2
         3
         4
         5
         6
         7
         8
         9
         10> for (i in 1:10) 
+ { 
+ resposta[i] = mean(round(rnorm(10, mean=7.5, sd=3.2)))- mean(round(rnorm(10, mean=6, sd=3))) 
+ } 
> resposta[i]
[1] 3.4
> 
> dif.total=simula(data.a,data.b,nsim=2000,teste="bi")
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 =  0.1 
> table(dif.total)
dif.total
  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 
 49 102 103  98  98  99 106  92  80  87  80  80  78  88  56  66  55  69  56  59 
  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 
 43  33  31  30  38  19  27  19  21  18  20  12  12  13  10   3   6   7   9   3 
  4 4.1 4.2 4.3 4.4 4.5 4.6 4.7 5.6 
  4   5   2   3   3   2   2   3   1 
> dif.depois=sum(dif.total>=3.4) 
> dif.depois
[1] 63
> p1=dif.depois/length(dif.total)
> p1
[1] 0.0315
> 

2-)Teste agora que a média da segunda amostra é maior que a 
primeira.

R:
resposta.2 = rep(NA, 10) 
for (j in 1:10) 
{ 
cat('\n\t', j) 
} 
for (j in 1:10) 
{ 
media.a=mean(round(rnorm(10, mean=7.5, sd=3.2)))
media.b=mean(round(rnorm(10, mean=6, sd=3)))
media.a
media.b
resposta.2[j] = media.a- media.b 
} 
resposta.2[j] 
dif.a.b = simula(data.a,data.b, nsim=2000, teste='uni') 
table(dif.a.b) 
#mostra o numero de ocorrencias de cada media 
menor.dif = sum(dif.a.b)<=dif.entre 
menor.dif 
dif.entre
p2 = menor.dif/length(dif.a.b) 
p2 


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

R:
t.test(data.a,data.b)

        Welch Two Sample t-test

data:  data.a and data.b 
t = -0.0612, df = 17.841, p-value = 0.9519
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -3.536558  3.336558 
sample estimates:
mean of x mean of y 
      6.1       6.2 


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

R:
par(mfrow=c(3,1))
hist(datas, freq=FALSE,xlim=c(1,20)) 
curve(exp=dnorm(x, mean=mean(data),sd=sd(data)),from=0,to=21, col="red", add=T) 
boxplot(datas~nomes.conj.) 
plot(1:20, datas, pch=rep(c(15,16), each=10), col=rep(1:2, each=10)) 
for (i in 1:10) 
{ 
lines(c(i,i),c(datas[i],media.a),col=1) 
} 
for(j in 11:20) 
{ 
lines(c(j,j),c(datas[j],media.b),col=2) 
} 
lines(c(1,10),c(media.a,media.a),col=1) 
lines(c(11,20),c(media.b,media.b),col=2) 
par(mfrow=c(1,1)) 

##Caixeta de NOVO?!


#Utilizando os dados da planilha caixeta.csv 

caixeta=read.table("caixeta.csv",header=T,as.is=T,sep=",")


1-)Calcule os valores de área basal por fuste 

R:
head(caixeta)
caixeta$parcela = as.factor(caixeta$parcela) 
caixeta$local = as.factor(caixeta$local) 
str(caixeta) 
caixeta$area.basal <- (caixeta$cap^2)/(4*pi)
caixeta$area.basal 
str(caixeta$area.basal) 
length(caixeta$fuste) 
head(caixeta)

2-)Calcule os valores de área basal por amostra em cada uma das 
localidades

R:
chauas = caixeta[(caixeta$local=="chauas"),] 
jureia = caixeta[(caixeta$local=="jureia"),] 
retiro = caixeta[(caixeta$local=="retiro"),] 
area.basal.chauas = tapply(chauas$area.basal, INDEX=chauas$parcela, FUN=sum) 
area.basal.jureia= tapply(jureia$area.basal, INDEX=jureia$parcela, FUN=sum) 
area.basal.retiro =tapply(retiro$area.basal, INDEX= retiro$parcela, FUN = sum) 
area.basal.chauas 
area.basal.jureia 
area.basal.retiro




3-)Produza gráficos para mostrar os dados

R:
par(cex.axis=1) 
par(cex.lab=1.3) 
par(bty="l")
boxplot(area.basal.chauas, area.basal.jureia, area.basal.retiro, ylab="Área basal", xlab="", main= "Área basal X Localidade") 
mtext(c("Chauas","Juréia","Retiro")
, side= 1 ,line=2,at=c(1,2,3)) 



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.

R:
#DESVIO QUADRÁTICO TOTAL
caixeta<-data.frame(basal.chauas,basal.jureia,basal.retiro)
caixeta
media.geral<-mean(c(basal.chauas,basal.jureia,basal.retiro))
media.geral
des.total<-caixeta-media.geral
des.total
des.total.q<-desvio.total^2
des.total.q
soma.des.total.q<-sum(des.total.q)
soma.des.total.q  

##DESVIO QUADRÁTICO INTRAGRUPOS
med.local<-mean(caixeta)
med.local
des.intra<-data.frame(basal.chauas,basal.jureia,basal.retiro)
des.intra$desvio.chauas<-caixeta$basal.chauas-med.local["basal.chauas"]
des.intra$desvio.jureia<-caixeta$basal.jureia-med.local["basal.jureia"]
des.intra$desvio.retiro<-caixeta$basal.retiro-med.local["basal.retiro"]
des.intra

desvio.chauas<-sum(desvios.intra["desvio.chauas"]^2)
desvio.chauas
desvio.jureia<-sum(desvios.intra["desvio.jureia"]^2)
desvio.jureia
desvio.retiro<-sum(desvios.intra["desvio.retiro"]^2)
desvio.retiro

soma.intra<-sum(desvio.chauas,desvio.jureia,desvio.retiro)
soma.intra   

##DESVIO QUADRÁTICO ENTRE GRUPOS
diff.entre<-soma.des.total.q-soma.intra
diff.entre    

## GRAUS DE LIBERDADE

grau.liberdade.total<-15-1   		
grau.liberdade.intra<-15-3  		
grau.liberdade.entre<-grau.liberdade.total-grau.liberdade.intra   
grau.liberdade.total
grau.liberdade.intra
grau.liberdade.entre

## DESVIO MÉDIO

desvio.medio.intra<-soma.intra/grau.liberdade.intra
desvio.medio.intra			
desvio.medio.entre<-soma.entre/grau.liberdade.entre
desvio.medio.entre
			

## CÁLCULO DE F

F<-desvio.medio.entre/desvio.medio.intra
F     

## CÁLCULO DE p

p<-pf(F, 2, 12, lower.tail = FALSE, log.p = FALSE)
p    


   FONTE      DESVIO QUADRÁTICO	GRAUS DE LIBERDADE	DESVIO MÉDIO	RAZÃO DAS VARIÂNCIAS	PROBABILIDADES
Entre grupos	1649495370    	      2    	         824747685	     1.329489	           0.3009299
Intra grupos	7444192546	     12	                 620349379		

Total	        9093687916           14			


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

R:
> area.basal<-c(basal.chauas,basal.jureia,basal.retiro) 
> locais <- as.factor(c(rep(c("basal.chauas","basal.jureia", "basal.retiro"), each=5))) 
> aov(area.basal~locais) 
Call:
   aov(formula = area.basal ~ locais)

Terms:
                    locais  Residuals
Sum of Squares  1649495370 7444192546
Deg. of Freedom          2         12

Residual standard error: 24906.81 
Estimated effects may be unbalanced


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

R:
> porc.=(diff.entre/soma.intra)*100
> porc.
[1] 22.15815
>