1 Uma senhora toma chá

Visite

#-----------------------------------------------------------------------
# O curioso caso da senhora que toma chá.

# Quantidade de maneiras de gerar dois grupos de 4 xícaras usando 8.
n_poss <- choose(8, 4)
n_poss

## [1] 70

# E se fossem 12 xícaras, 6 de cada?
choose(12, 6)

## [1] 924

# De 4 xícaras selecionadas para um grupo, X representa o número de
# acertos.
n_acertar <- sapply(4:0,
                    FUN = function(i) {
                        # (formas de acertar x em 4) * (formas de errar x em 4).
                        choose(4, i) * choose(4, 4 - i)
                    })
n_acertar

## [1]  1 16 36 16  1

# Pr(Acerto total) = Pr(X = 4) = 1/70 < 5%.
cbind(X = 4:0, "Pr(X = x)" = n_acertar/n_poss)

##      X  Pr(X = x)
## [1,] 4 0.01428571
## [2,] 3 0.22857143
## [3,] 2 0.51428571
## [4,] 1 0.22857143
## [5,] 0 0.01428571

# Comprovando por simulação.
# Uma particular configuração das xícaras.
v <- rep(0:1, each = 4)

# Com qual frequência, ao acaso, a mesma configuração ocorre?
mean(replicate(n = 100000,
               expr = all(sample(v) == v)))

## [1] 0.01427

# Essa forma de fazer é mais realista mas o resultado é o memso.
mean(replicate(n = 100000,
               expr = all(sample(v) == sample(v))))

## [1] 0.01398

2 Teste para a diferença de médias

#-----------------------------------------------------------------------
# Exemplo com teste para a diferença de médias.

# Comprimentos de mandíbulas de cães pré-históricos.
m <- c(120, 107, 110, 116, 114, 111, 113, 117, 114, 112) # machos.
f <- c(110, 111, 107, 108, 110, 105, 107, 106, 111, 111) # fêmeas.

#  Distribuição acumulada empírica.
plot(ecdf(m), xlim = range(c(m, f)), col = "cyan")
lines(ecdf(f), col = "magenta")
rug(m, col = "cyan")
rug(f, col = "magenta")

# Diferença de média observada.
d <- mean(m) - mean(f)
d

## [1] 4.8

# Todos as combinações possíveis entre 20 objetos tomados 10 a 10.
choose(n = 20, k = 10)

## [1] 184756

#-----------------------------------------------------------------------
# Solução com todas as permutações possíveis (exaustão).

# Para construir todas as combinações possíveis.
k <- combn(x = 1:20, m = 10)
dim(k)

## [1]     10 184756

# Vetor com os valores dos dois grupos concatenados.
mf <- c(m, f)

# Vetor de zeros.
g <- integer(20)

# Calcula a diferença para cada uma das permutações.
D <- apply(k,
           MARGIN = 2,
           FUN = function(i) {
               # Alguns zeros são convertidos para 1 criando 2 grupos.
               g[i] <- 1L
               # Cálculo da estatística de teste.
               -diff(tapply(mf, g, FUN = mean))
           })

# Histograma da distribuição da estatística sob H_0.
hist(D, col = "gray50")
rug(D)
abline(v = d, col = 2)

# Distribuição acumulada empírica da estatística sob H_0.
plot(ecdf(D), cex = 0)
rug(D)
abline(v = d, col = 2)

# P-valor do teste (bilateral).
2 * sum(D >= d)/length(D)

## [1] 0.003334127

#-----------------------------------------------------------------------
# Com reamostragem sem reposição (não necessáriamente exaustivo).

# Variáveis que indentifica os grupos.
g <- rep(1:2, each = 10)

# Apenas para conferir.
cbind(g, mf)

##       g  mf
##  [1,] 1 120
##  [2,] 1 107
##  [3,] 1 110
##  [4,] 1 116
##  [5,] 1 114
##  [6,] 1 111
##  [7,] 1 113
##  [8,] 1 117
##  [9,] 1 114
## [10,] 1 112
## [11,] 2 110
## [12,] 2 111
## [13,] 2 107
## [14,] 2 108
## [15,] 2 110
## [16,] 2 105
## [17,] 2 107
## [18,] 2 106
## [19,] 2 111
## [20,] 2 111

# Replicando a diferença para grupos formados por aleatorização.
D <- replicate(9999, {
    gg <- sample(g)
    -diff(tapply(mf, gg, FUN = mean))
})

# Tem que juntar a estatística observada com as simuladas.
D <- c(D, d)

hist(D, col = "gray50")
rug(D)
abline(v = d, col = 2)

plot(ecdf(D), cex = 0)
rug(D)
abline(v = d, col = 2)

# P-valor do teste (bilateral).
2 * sum(D >= d)/length(D)

## [1] 0.0054

3 Teste para a correlação

#-----------------------------------------------------------------------
# Teste de aleatorização para a correlação.

# N = 5 para um par de medidas.
x <- c(4.1, 8.3, 2.9, 10.8, 9.5)
y <- c(3.7, 5.1, 1.0, 7.7, 8.9)
cbind(x, y)

##         x   y
## [1,]  4.1 3.7
## [2,]  8.3 5.1
## [3,]  2.9 1.0
## [4,] 10.8 7.7
## [5,]  9.5 8.9

plot(x, y)

# Estatística de teste na amostra original: coeficiente de correlação de
# Pearson.
r0 <- cor(x, y, method = "pearson")
r0

## [1] 0.9228669

# Todas as permutações possiveis: 5! = 120 do vetor x.
X <- gtools::permutations(n = length(x), r = length(x), v = x)
str(X)

##  num [1:120, 1:5] 2.9 2.9 2.9 2.9 2.9 2.9 2.9 2.9 2.9 2.9 ...

head(X)

##      [,1] [,2] [,3] [,4] [,5]
## [1,]  2.9  4.1  8.3  9.5 10.8
## [2,]  2.9  4.1  8.3 10.8  9.5
## [3,]  2.9  4.1  9.5  8.3 10.8
## [4,]  2.9  4.1  9.5 10.8  8.3
## [5,]  2.9  4.1 10.8  8.3  9.5
## [6,]  2.9  4.1 10.8  9.5  8.3

# As 120 correlações obtidas para cada arranjo.
r <- apply(X, MARGIN = 1, FUN = cor, y = y, method = "spearman")

#  Gráfico da distribuição da estatística de teste.
plot(ecdf(r), cex = 0)
rug(r)
abline(v = r0, col = 2)

# P-valor do teste.
2 * sum(r >= r0)/length(r)

## [1] 0.01666667

4 Índice de Moran

#-----------------------------------------------------------------------
# Índice de Moran para medir dependência espacial.

# Coordenadas dos eventos em uma malha regular 8 x 8.
x <- 1:8
y <- 1:8

# Construção da matriz de pesos que determina a vizinhança entre
# observações.
ind <- expand.grid(i = 1:length(x),
                   j = 1:length(y))
#  Função que determina o peso entre duas localizações na malha.
f <- function(i, j) {
    u <- min(3, sum(abs(ind[i, ] - ind[j, ])))
    w <- c(0, 1, sqrt(1/2), 0)[u + 1]
    return(w)
}
#  Cria os pesos, matriz (8^2) x (8^2) = 64 x 64.
w <- matrix(0, nrow = nrow(ind), ncol = nrow(ind))
for (i in 1:nrow(ind)) {
    for (j in 1:nrow(ind)) {
        w[i, j] <- f(i, j)
    }
}
#  Normaliza.
w <- w/sum(w)

#  Gráfico. Valores claros indicam maior peso entre observações.
image(w, asp = 1, col = gray.colors(100))

#  Lógica do índica de Moran: correlação entre valores observados e
#  média dos vizinhos. Exemplo com valores simulados.
xx <- rnorm(64)
cor(cbind("Valores observados" = xx,
          "Média dos vizinhos" = as.vector(xx %*% w)))

##                    Valores observados Média dos vizinhos
## Valores observados          1.0000000         -0.1279173
## Média dos vizinhos         -0.1279173          1.0000000

# Índice de Moran.
moran <- function(x, weights) {
    # Tamanho da amostra.
    n <- length(x)
    # Valores normalizados.
    z <- as.vector((x - mean(x))/sd(x))
    # Índice de Moran.
    as.vector(z %*% weights %*% (z * sqrt(n/(n - 1))))
}

# Teste de permutação com saída gráfica.
ppt <- function(z, w, N = 10000, stat, ...) {
    # Índice de Moran por reamostragem.
    sim <- replicate(N,
                     moran(sample(z), w))
    # Determina o p-valor.
    p.value <- mean((all <- c(stat, sim)) >= stat)
    # Histograma da distribuição empírica sob H_0.
    hist(sim,
         sub = paste("p =", round(p.value, 4)),
         xlim = range(all),
         ...)
    abline(v = stat, col = "#903030", lty = 3, lwd = 2)
    return(p.value)
}

# Observações simuladas.
set.seed(17)
par(mfrow = c(2, 3))

# Dados com dependência espacial ---------------------------------------
# Indução de autocorrelação por meio de uma tendência.
z <- matrix(rexp(length(x) * length(y),
                 outer(x, y^2)),
            length(x))
image(log(z), main = "Com dependência")

cor(cbind("Valores observados" = as.vector(z),
          "Média dos vizinhos" = as.vector(as.vector(z) %*% w)))

##                    Valores observados Média dos vizinhos
## Valores observados          1.0000000          0.1335676
## Média dos vizinhos          0.1335676          1.0000000

# Índice de Moran com dados originais.
stat <- moran(z, w)
stat

## [1] 0.06499871

hist(z)
ppt(z, w, stat = stat, main = "I de Moran", xlab = "I")

## [1] 0.01459854

# Dados sem dependência espacial ---------------------------------------
# Geração de de um conjunto de dados sob hipótese nula.
z <- matrix(rnorm(length(x) * length(y), 0, 1/2), length(x))
image(z, main = "Sem dependência")

cor(cbind("Valores observados" = as.vector(z),
          "Média dos vizinhos" = as.vector(as.vector(z) %*% w)))

##                    Valores observados Média dos vizinhos
## Valores observados         1.00000000        -0.04639892
## Média dos vizinhos        -0.04639892         1.00000000

# Índice de Moran com dados originais.
stat <- moran(z, w)
stat

## [1] -0.0121901

hist(z)
ppt(z, w, stat = stat, main = "I de Moran", xlab = "I")

## [1] 0.4478552

Estatística Computacional II

`leg.ufpr.br/~walmes/ensino/EC2`

Aplicações de testes de aleatorização

Prof. Walmes M. Zeviani

2018-09-10

1 Uma senhora toma chá

2 Teste para a diferença de médias

3 Teste para a correlação

4 Índice de Moran