Técnicas de Monte Carlo

Introdução

Simulação e aleatorizações baseadas em dados ou distribuições teóricas buscando soluções numéricas.

• teste de hipóteses
• medidas de precisão de estimativas
• otimizadores
• integração numérica
• algoritmos de amostragem

Definições

Monte Carlo

Técnicas de simulação buscando resultado numérico

• simulações aleatórias (numérico)
• distribuição conhecida (normal, poisson)
• MCMC “Markov Chain Monte Carlo”

Reamostragem

Técnicas de reamostragem de dados

• precisão de estimativa (bootstrap, jackknife)
• teste de significância (reordenação, permutação)
• validação de modelos (subconjuntos)

Definições

Teste de Permutação

• reordenamento (rótulos) em todas as combinações possíveis (teste exato de Fisher)
• combinações possíveis de 10 valores:

factorial(10)
[1] 3628800

• Teste de Monte Carlo: uma amostra das combinações

Reposição

Bootstrap

Reamostragem com reposição

Jackknife

Reamostragem de subconjunto

Inferência: precisão de uma estimativa

Implicações

Não assumem a distribuição de probabilidade teórica

Vantagens

• poucas restrições (dados)
• intuitiva (conhecimento matemático)
• poucos pressupostos
• assintótica

Desvantagens

• dificuldades computacionais
• resultado pode variar (estatítica de interesse)
• domínio de inferência restrito
• necessita cenário nulo adequado (complexo)

Teste de Hipóteses

1. Definir a estatística de interesse (EI)
2. Estabelecer o cenário nulo
3. Reamostrar, reordenar ou simular o cenário nulo
4. Calcular a EI no cenário nulo
5. Criar a distribuição dos pseudovalores da EI
6. Posicionar o observado na distribuição dos pseudovalores
7. Calcular o p-valor

Euterpe edulis

amospalm <- read.table("amostraEut.txt", header = TRUE, as.is=TRUE, sep = "\t")
str(amospalm)

'data.frame':   200 obs. of  2 variables:
 $ local: chr  "PECB" "PECB" "PECB" "PECB" ...
 $ pap  : num  46 29 48 42 21 25 49 41 26 47 ...

local <- factor(amospalm$local, levels = c("PEIC", "PECB"))
tapply(amospalm$pap, local, mean)

  PEIC   PECB 
25.495 35.640 

(difpalm <- diff(tapply(amospalm$pap, local,mean)))

  PECB 
10.145 

Representação dos dados

par(mar= c(4,4,2,2),las=1, cex =1.5,cex.axis = 1.2, cex.lab= 1.2)
boxplot(amospalm$pap ~ local,  col="gray", xlab = "Local",
        ylab = "PAP (cm)", outline=FALSE) 
points(x= jitter(rep(1,100), factor=2), y = amospalm$pap[local=="PEIC"],
       pch=16, col=rgb(1,0,0, alpha=0.3) )  
points(x= jitter(rep(2,100), factor=2), y = amospalm$pap[local=="PECB"],
       pch=16, col=rgb(0,0,1, alpha=0.3)) 
mtext(at=1.5, side=3, line=0.5, text =paste("diferença =",
                                            round(difpalm,1)), cex = 2)

Representação dos dados

Cenário Nulo

A diferença observada pode ser gerada pelo acaso?

Cenário Nulo

> cenulo <- sample(amospalm$pap)
> tapply(cenulo, local,mean)

  PEIC   PECB 
31.465 29.670 

> diff(tapply(cenulo,local,mean))

  PECB 
-1.795 

> diff(tapply(sample(amospalm$pap), local, mean))

 PECB 
2.425 

> diff(tapply(sample(amospalm$pap), local, mean))

  PECB 
-2.245 

Distribuição Nula

> dnulo <- rep(NA, 1000)
> dnulo[1] <- diff(tapply(amospalm$pap, local, mean))
> for(i in 2:1000)
+ {
+ dnulo[i] <- diff(tapply(sample(amospalm$pap), local, mean))
+ }

Null Anima

Distribuição Nula

Incerteza Sobre a Afirmação

A probabilidade da diferença observada ter sido gerada pelo acaso.

> sum(dnulo >= difpalm) 

[1] 1

> sum(dnulo >= difpalm)/length(dnulo)

[1] 0.001

Incerteza Sobre a Afirmação

A probabilidade da diferença observada ter sido gerada pelo acaso.

> sum(dnulo >= difpalm) 

[1] 1

> sum(dnulo >= difpalm)/length(dnulo)

[1] 0.001

p-valor

Distribuição espacial

Plantas estão distribuídas aleatoriamente em uma parcela?

Estatística de interesse (EI)

• média da distância ao vizinho mais próximo

Estatística de interesse (EI)

• média da distância ao vizinho mais próximo

xy[1:3,]

  xp yp
1  8  3
2  6 14
3  5 29

distmat <- as.matrix(dist(xy, diag = FALSE, upper=TRUE))

Estatística de interesse (EI)

distmat[1:3, 1:3]

         1        2       3
1  0.00000 11.18034 26.1725
2 11.18034  0.00000 15.0333
3 26.17250 15.03330  0.0000

diag(distmat) <- NA

Cálculo da estatística de interesse

diag(distmat) <- NA
nnd<- apply(distmat, 1, min, na.rm=TRUE)
tail(nnd)

      95       96       97       98       99      100 
2.000000 6.324555 6.708204 2.828427 5.099020 4.000000 

mean(nnd)

[1] 4.403911

Definir cenário nulo

• completa aleatoriedade espacial

Simular o cenário nulo

rx <- round(runif(100,0,100),0)
ry <- round(runif(100,0,100),0)
rxy <- data.frame(rx, ry)

Calcular a EI no nulo

rdistmat <- as.matrix(dist(rxy, diag = FALSE, upper=TRUE)) 
diag(rdistmat) <- NA
rnnd <- apply(rdistmat, 1, min, na.rm=TRUE)
mean(rnnd)

[1] 5.938513

Calcular a EI no nulo

• completa aleatoriedade espacial

Distribuição da EI no nulo

• definir o número de simulações
• criar o objeto de resultado das simulações

nsim = 1000
cnulo = rep(NA, nsim)
cnulo[1] <- mean(nnd)
cnulo[1:5]

[1] 4.403911       NA       NA       NA       NA

Distribuição da EI no nulo

• criar o ciclo
• armazenar o resultado na posição

for(i in 2:nsim)
{
rx <- round(runif(100,0,100),0)
ry <- round(runif(100,0,100),0)
rxy <- data.frame(rx, ry)
rdistmat <- as.matrix(dist(rxy, diag = FALSE, upper=TRUE)) 
diag(rdistmat) <- NA
rnnd <- apply(rdistmat, 1, min, na.rm=TRUE)
cnulo[i] <- mean(rnnd)
}

Simulação

Calcular o p-valor

cnulo[1:10]

 [1] 4.403911 5.356657 5.147929 5.214187 4.814227 5.571546 5.211867
 [8] 5.604144 5.046847 5.423152

ccnulo <- cnulo - mean(cnulo)

Calcular o p-valor

ccnulo <- cnulo - mean(cnulo)
hist(ccnulo)
abline(v= ccnulo[1], lty=2, col="red")

Calcular o p-valor

(nmaior <- sum(abs(ccnulo) >= abs( ccnulo[1])))

[1] 7

(pval <- nmaior/length(ccnulo))

[1] 0.007

Resultado:

• pontos mais próximos do que o esperado pelo cenário de completa aleatoriedade espacial

Completa Aleatoriedade Espacial

Categórica com 3 níveis

are=c(6,10,8,6,14,17,9,11,7,11)
arg=c(17,15,3,11,14,12,12,8,10,13)
hum=c(13,16,9,12,15,16,17,13,18,14)
crop <- data.frame(solo = rep(c("are", "arg", "hum"), each=10), colhe = c(are, arg, hum))
head(crop)

  solo colhe
1  are     6
2  are    10
3  are     8
4  are     6
5  are    14
6  are    17

Anova: partição da variação

\[F=\frac{\sigma_{entre}^2}{\sigma_{intra}^2}\]

Estatística de interesse

Médias dos solos

(msolo <- tapply(crop$colhe, crop$solo, mean))

 are  arg  hum 
 9.9 11.5 14.3 

Média Geral

(mcolhe <- mean(crop$colhe))

[1] 11.9

Estatística de interesse

Soma das diferenças

msolo - mcolhe

 are  arg  hum 
-2.0 -0.4  2.4 

sum(msolo - mcolhe)

[1] 0

(difobs <- sum(abs(msolo-mcolhe)))

[1] 4.8

Cenário Nulo

Distribuição da EI no nulo

• definir o número de simulações
• criar o objeto de resultado das simulações

nsim = 1000
difnulo = rep(NA, nsim)
difnulo[1] <- difobs
difnulo[1:5]

[1] 4.8  NA  NA  NA  NA

Distribuição da EI no nulo

• criar o ciclo
• armazenar o resultado na posição

for(i in 2:nsim)
{
scolhe <- sample(crop$colhe)
smsolo <- tapply(scolhe, crop$solo, mean)
difnulo[i] <- sum(abs(smsolo - mcolhe))
}

Simulação

Calcular o p-valor

difnulo[1:10]

 [1] 4.8 0.8 1.6 3.0 3.0 5.4 3.0 0.4 2.4 3.4

sum(difnulo >= difnulo[1])

[1] 27

sum(difnulo >= difnulo[1])/length(difnulo)

[1] 0.027

Anova Curso R

Davis (1990). Appetite (15)13-21

library(car)
data(Davis)
summary(Davis)

 sex         weight          height          repwt            repht      
 F:112   Min.   : 39.0   Min.   : 57.0   Min.   : 41.00   Min.   :148.0  
 M: 88   1st Qu.: 55.0   1st Qu.:164.0   1st Qu.: 55.00   1st Qu.:160.5  
         Median : 63.0   Median :169.5   Median : 63.00   Median :168.0  
         Mean   : 65.8   Mean   :170.0   Mean   : 65.62   Mean   :168.5  
         3rd Qu.: 74.0   3rd Qu.:177.2   3rd Qu.: 73.50   3rd Qu.:175.0  
         Max.   :166.0   Max.   :197.0   Max.   :124.00   Max.   :200.0  
                                         NA's   :17       NA's   :17     

Davis (1990). Appetite (15)13-21

Davis[Davis$weight >150,]

   sex weight height repwt repht
12   F    166     57    56   163

Davis <- Davis[-12, 1:3]
str(Davis)

'data.frame':   199 obs. of  3 variables:
 $ sex   : Factor w/ 2 levels "F","M": 2 1 1 2 1 2 2 2 2 2 ...
 $ weight: int  77 58 53 68 59 76 76 69 71 65 ...
 $ height: int  182 161 161 177 157 170 167 186 178 171 ...

lmdavis<- lm(weight ~ height, data =Davis)
coef(lmdavis)

(Intercept)      height 
-130.746984    1.149222 

REGRESSÃO: gráfico

REGRESSÃO: estatística de interesse

• inclinação da reta (\(\beta\))

\[ \hat{\beta} =\frac{ \sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})}{ \sum_{i=1}^{n}(x_i - \bar{x})^2}\]

REGRESSÃO: cenário nulo

Davis$height[1:5]

[1] 182 161 161 177 157

sample(Davis$height)[1:5]

[1] 175 173 157 175 179

Davis$simh <-  sample(Davis$height)
kable(head(Davis))

Permutando a altura

coef(lmdavis)

(Intercept)      height 
-130.746984    1.149222 

bobs <- round(coef(lmdavis)[2],3)

lmsim <- lm(weight ~ simh, data = Davis) 
(bsim <- round(coef(lmsim)[2], 3))

  simh 
-0.136 

Cenário Nulo

slope simulado = -0.136

Distribuição de pseudovalores

nSim = 1000
psv = rep(NA, nSim )
psv[1] <- coef(lm(weight ~ height, data=Davis))[[2]]
for(i in 2:nSim)
{
psv[i] <- coef(lm(weight ~ sample(height), data=Davis))[[2]]
}

Simulação

REGRESSÃO: p-valor

pvalor = sum(abs(psv) >= abs(psv[1]))/length(psv)
pvalor

[1] 0.001

Regressão

ANCOVA por reamostragem

Que perguntas?

ANCOVA por reamostragem

Que perguntas podemos fazer?

1. há relação entre peso e altura
2. a relação entre os sexos é a mesma, mas há um efeito de ser macho?
3. os sexos apresentam relações diferentes?

ANCOVA por reamostragem

a relação é a mesma, mas há um efeito de ser macho:

\[ \hat{\alpha}_m \neq \hat{\alpha}_f\]

ANCOVA por reamostragem

Os sexos apresentam relações diferentes:

\[ \hat{\beta}_m \neq \hat{\beta}_f\]

sample(.., replace = TRUE)

ae <- letters[1:5]
sample(ae)

[1] "a" "e" "d" "c" "b"

sample(ae, replace=TRUE)

[1] "d" "b" "b" "a" "d"

sample(ae, size=10, replace=TRUE)

 [1] "a" "e" "d" "b" "b" "c" "a" "a" "c" "d"

Bootstrap

Intervalo de confiança por percentil

macho=c(120,107,110,116, 114, 111, 113,117,114,112)
(media.m=mean(macho))

[1] 113.4

Chacal macho: resultados

Qual a minha confiança sobre uma estimativa?

intervalo de confiança da média

macho

 [1] 120 107 110 116 114 111 113 117 114 112

sample(macho, replace = TRUE)

 [1] 114 116 107 117 110 112 120 112 117 111

sample(macho, replace = TRUE)

 [1] 112 117 116 116 111 107 114 116 107 117

Estimativa bootstrap

mean(macho)

[1] 113.4

mean(sample(macho, replace = TRUE))

[1] 112.2

mean(sample(macho, replace = TRUE))

[1] 111.3

mean(sample(macho, replace = TRUE))

[1] 115.6

Bootstrap

nBoot = 1000
bvalue <- rep(NA, nBoot )
bvalue[1] <- mean(macho)
  for(i in 2:nBoot)
  {
    bvalue[i] <- mean(sample(macho, replace = TRUE))
  }

Bootstrap

Chacal macho: Resultado

Rsampling

https://github.com/lageIBUSP/Rsampling-shiny

install.packages(c("Rsampling", "shiny", "PerformanceAnalytics"))
library("shiny")
runApp("/home/aao/Ale2016/AleCursos/Planejamento&Analise/Rsampling-shiny-1.6.1")

Bibliografia

Atividades da Tarde