Adiciona scripts 4, 5 e 6.

parent 1dd11794
#-----------------------------------------------------------------------
# Teste de aleatórização.
# Comprimentos de crânios de cães pré históricos.
m <- c(120, 107, 110, 116, 114, 111, 113, 117, 114, 112)
f <- c(110, 111, 107, 108, 110, 105, 107, 106, 111, 111)
# Diferença de média.
d <- mean(m) - mean(f)
d
# Aplicando um teste t.
t.test(x = m, y = f, var.equal = TRUE)
# Vetor com os valores dos dois grupos.
mf <- c(m, f)
# Variáveis que indentifica os grupos.
g <- rep(1:2, each = 10)
cbind(mf, g)
# Replicando a diferença para grupos formados por aleatorização.
D <- replicate(999, {
y <- sample(mf, size = length(mf), replace = FALSE)
-diff(tapply(y, g, FUN = mean))
})
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.
2 * sum(D >= d)/length(D)
#-----------------------------------------------------------------------
# Teste para 3 grupos.
# Valores observados (y) e grupos (g).
y <- c(7, 8, 11, 12, 14, 19, 21, 122, 0, 2, 5, 9)
g <- gl(3, 4)
cbind(y, g)
plot(y ~ as.numeric(g))
# Valor observado da estatística F.
fobs <- anova(lm(y ~ g))[1, "F value"]
fsim <- replicate(9999, {
y <- sample(y, size = 12, replace = FALSE)
anova(lm(y ~ g))[1, "F value"]
})
fsim <- c(fsim, fobs)
# P-valor do teste.
sum(fsim >= fobs)/length(fsim)
#-----------------------------------------------------------------------
# Teste de aleatorização para a correlação.
x <- c(4, 8, 2, 10, 9)
y <- c(3, 5, 1, 7, 8)
r0 <- cor(x, y)
r0 <- cor(x, y, method = "spearman")
library(gtools)
# Todas as permutações possiveis = 5! = 120.
Y <- permutations(n = length(y), r = length(y), v = y)
str(Y)
r <- apply(Y, MARGIN = 1, FUN = cor, y = x)
# P-valor do teste.
sum(r >= r0)/length(r)
#-----------------------------------------------------------------------
# Índice de Moran para medir dependência espacial.
# http://gis.stackexchange.com/questions/161887/significance-test-for-morans-i-using-monte-carlo-simulation
browseURL("https://en.wikipedia.org/wiki/Moran's_I")
# Moran's I
# https://en.wikipedia.org/wiki/Moran's_I
# Coordenadas dos eventos.
x <- 1:8
y <- 1:8
# Matriz de pesos.
ind <- expand.grid(i = 1:length(x), j = 1:length(y))
f <- function(i, j) {
u <- min(3, sum(abs(ind[i, ] - ind[j, ])))
c(0, 1, sqrt(1/2), 0)[u + 1]
}
w <- matrix(0, nrow(ind), nrow(ind))
for (i in 1:nrow(ind)) {
for (j in 1:nrow(ind)) {
w[i, j] <- f(i, j)
}
}
w <- w/sum(w)
image(w)
# Índice de Moran.
moran <- function(x, weights) {
n <- length(x)
z <- as.vector((x - mean(x))/sd(x))
as.vector(z %*% weights %*% (z * sqrt(n/(n - 1))))
}
# Teste de permutação com gráfico.
ppt <- function(z, w, N = 10000, ...) {
stat <- moran(z, w)
sim <- replicate(N, moran(sample(z, length(z)), w))
p.value <- mean((all <- c(stat, sim)) >= stat)
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))
# 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 = "Autocorrelated Data")
hist(z)
ppt(z, w, main = "Null Distribution of Moran's I", xlab = "I")
# 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 = "Uncorrelated Data")
hist(z)
ppt(z, w, main = "Null Distribution of Moran's I", xlab = "I")
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
# Aplicação de Jackknife em modelos de regressão não linear.
library(latticeExtra)
library(car)
library(alr3)
str(segreg)
help(segreg, help_type = "html")
xyplot(C ~ Temp, data = segreg, type = c("p", "smooth"))
# Curva "palpite".
start <- list(th0 = 75, th1 = 0.5, th2 = 50)
xyplot(C ~ Temp, data = segreg) +
layer(panel.curve(th0 + th1 * (x - th2) * (x >= th2) +
0 * (x < th2), lty = 2),
data = start)
# Ajuste.
n0 <- nls(C ~ th0 + th1 * (Temp - th2) * (Temp >= th2) +
0 * (Temp < th2),
data = segreg,
start = start)
# Estimativas e medidas de ajuste.
summary(n0)
# Intervalos de confiança.
# confint(n2)
confint.default(n0)
# Observados e preditos.
xyplot(C ~ Temp, data = segreg) +
layer(panel.curve(th0 + th1 * (x - th2) * (x >= th2) +
0 * (x < th2)),
data = as.list(coef(n0)))
#-----------------------------------------------------------------------
# Estimativas Jackknife para os parâmetros.
n <- nrow(segreg)
theta <- coef(n0)
thetaj <- matrix(0, nrow = n, ncol = length(theta))
colnames(thetaj) <- names(theta)
# Ajustes deixando uma observação de fora.
for (i in 1:n) {
n1 <- nls(C ~ th0 + th1 * (Temp - th2) * (Temp >= th2) +
0 * (Temp < th2),
data = segreg[-i, ],
start = coef(n0))
thetaj[i, ] <- n * theta - (n - 1) * coef(n1)
}
thetaj
pairs(thetaj)
# Estimativas pontuais.
cbind(Jack = colMeans(thetaj),
MLE = theta)
# Erros padrões.
cbind(Jack = apply(thetaj, MARGIN = 2, sd)/sqrt(n),
MLE = summary(n0)$coefficients[, "Std. Error"])
#-----------------------------------------------------------------------
# Aplicação na Mediana (função "discretas").
x <- sort(precip)
M <- median(x)
M
n <- length(x)
i <- 1:length(x)
y <- sapply(i,
FUN = function(i) {
ep <- median(x[-i])
pv <- n * M - (n - 1) * ep
return(pv)
})
stem(y)
# NOTE: Explique porque o Jackknife para a media só retorna 3 valores?
#-----------------------------------------------------------------------
# Capacidade de predição.
str(cars)
plot(dist ~ speed, cars)
i <- 1:nrow(cars)
# As funções interno() e externo() retornam as correlações entre predito
# e observado tendo a observação predita dentro da amostra (interno) ou
# fora (leave-one-out).
externo <- function(grau) {
Y <- sapply(i,
FUN = function(i) {
m0 <- lm(dist ~ poly(speed, grau),
data = cars[-i, ])
yi <- predict(m0, newdata = cars[i, ])
return(cbind(yi, cars[i, ]$dist))
})
# R² do observado com o predito leave-one-out.
cor(Y[1, ], Y[2, ])^2
}
interno <- function(grau) {
# Ajuste do modelo com todos os dados.
m0 <- lm(dist ~ poly(speed, grau), data = cars)
summary(m0)$r.squared
}
g <- 1:5
int <- sapply(g, interno)
ext <- sapply(g, externo)
xyplot(int + ext ~ g,
type = "o",
auto.key = TRUE)
#-----------------------------------------------------------------------
# Inferência para a DL 50.
library(labestData)
data(PaulaTb3.12)
str(PaulaTb3.12)
help(PaulaTb3.12, h = "html")
# Converte a resposta para binário.
PaulaTb3.12$resp <- as.integer(PaulaTb3.12$resp) - 1
plot(resp ~ vol, data = PaulaTb3.12)
m0 <- glm(resp ~ vol,
data = PaulaTb3.12,
family = binomial)
summary(m0)
dl <- -coef(m0)[1]/coef(m0)[2]
# str(m0$family)$link
plot(resp ~ vol, data = PaulaTb3.12)
curve(m0$family$linkinv(coef(m0)[1] + coef(m0)[2] * x),
add = TRUE)
abline(v = dl, h = 0.5, lty = 2)
n <- nrow(PaulaTb3.12)
i <- 1:n
DL <- sapply(i,
FUN = function(i) {
m0 <- glm(resp ~ vol,
data = PaulaTb3.12[-i, ],
family = binomial)
# Estimativa Parcial.
DL <- -coef(m0)[1]/coef(m0)[2]
# Pseudo-valor.
n * dl - (n - 1) * DL
})
m <- mean(DL)
e <- sd(DL)/sqrt(n)
plot(density(DL))
rug(DL)
#-----------------------------------------------------------------------
# Tempo do primeiro lugar nas provas de Monza desde 1950.
# Lendo tabela direto do site.
da <- XML::readHTMLTable("http://www.4mula1.ro/s/?cdx=9b0870e&pos=1&idt=44")
da <- da[[1]]
summary(da)
# Anos das corridas.
da$ano <- as.integer(sub("^(\\d{4})-.*$", "\\1", da$Date))
# Tempo do primeiro colocado na prova.
library(lubridate)
x <- hms(da$Time)
da$tpro <- 60 * hour(x) + minute(x) + second(x)/60
#-----------------------------------------------------------------------
# browseURL("http://www.4mula1.ro/statistics/grandprix/")
xyplot(tpro ~ ano, data = da, type = c("p", "smooth"))
# Cria uma variável que começa em 0.
da$anoc <- da$ano - min(da$ano)
da <- na.omit(da)
plot(tpro ~ anoc, data = da)
start <- list(b0 = 160, b1 = -3, xb = 25)
with(start, {
curve(b0 + b1 * x * (x <= xb) + b1 * xb * (x > xb),
add = TRUE, col = 2)
})
# Ajuste do modelo não linear.
n0 <- nls(tpro ~ b0 +
b1 * anoc * (anoc <= xb) +
xb * b1 * (anoc > xb),
data = da, start = start)
# Resumo do ajuste.
summary(n0)
# Intervalos de confiança.
# confint(n0)
confint.default(n0)
# Curva ajustada.
plot(tpro ~ anoc, data = da)
with(as.list(coef(n0)), {
curve(b0 + b1 * x * (x <= xb) + b1 * xb * (x > xb),
add = TRUE,
col = 4)
})
#-----------------------------------------------------------------------
# Estimativas Jackknife para os parâmetros.
n <- nrow(da)
theta <- coef(n0)
thetaj <- matrix(0, nrow = n, ncol = length(theta))
colnames(thetaj) <- names(theta)
for (i in 1:n) {
n1 <- nls(tpro ~ b0 +
b1 * anoc * (anoc <= xb) +
xb * b1 * (anoc > xb),
data = da[-i, ],
start = coef(n0))
thetaj[i, ] <- n * theta - (n - 1) * coef(n1)
}
thetaj
# pairs(M)
pairs(thetaj)
jktheta <- colMeans(thetaj)
cbind(theta, jktheta)
jkse <- apply(thetaj, MARGIN = 2, sd)/sqrt(n)
cbind(se = summary(n0)$coefficients[, 2], jkse = jkse)
jktheta + pnorm(0.975) * jkse
jkic <- sweep(outer(jkse,
c(fit = 0, lwr = -1, upr = 1) * pnorm(0.975),
FUN = "*"),
MARGIN = 1, STATS = jktheta, FUN = "+")
jkic
#-----------------------------------------------------------------------
library(nlstools)
ls("package:nlstools")
j0 <- nlsJack(n0)
summary(j0)
help(nlsJack, h = "html")
# Ainda não sei reproduzir o IC calculado pela nlsJack mas está descrito
# em Seber & Wild (1989).
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
library(alr3)
str(turk0)
plot(Gain ~ A, data = turk0)
# Lista com os valores iniciais.
start <- list(Int = 625, Ass = 180, Mei = 0.1)
# Adiciona as curvas.
with(start, {
curve(Int + Ass * A/(Mei + A),
xname = "A", add = TRUE, col = 2)
curve(Int + Ass * (1 - 2^(-A/Mei)),
xname = "A", add = TRUE, col = 4)})
# Ajuste do primeiro modelo.
n0 <- nls(Gain ~ Int + Ass * A/(Mei + A),
data = turk0,
start = start)
# Ajuste do segundo modelo.
n1 <- nls(Gain ~ Int + Ass * (1 - 2^(-A/Mei)),
data = turk0,
start = start)
cbind(coef(n0), coef(n1))
#-----------------------------------------------------------------------
n <- nrow(turk0)
i <- 1:n
myboot0 <- function() {
n0 <- nls(Gain ~ Int + Ass * A/(Mei + A),
data = turk0[sample(i, n, replace = TRUE), ],
start = start)
coef(n0)
}
myboot1<- function() {
n0 <- nls(Gain ~ Int + Ass * (1 - 2^(-A/Mei)),
data = turk0[sample(i, n, replace = TRUE), ],
start = start)
coef(n0)
}
#-----------------------------------------------------------------------
b0 <- replicate(999, myboot0())
b1 <- replicate(999, myboot1())
b0 <- cbind(b0, coef(n0))
b1 <- cbind(b1, coef(n1))
str(b0)
# Vício.
cbind(n0 = rowMeans(b0) - coef(n0),
n1 = rowMeans(b1) - coef(n1))
# Variância.
cbind(n0 = apply(b0, MARGIN = 1, var),
n1 = apply(b1, MARGIN = 1, var))
# Função para calcular o erro quadrático médio.
eqm <- function(x) {
sum((x - tail(x, 1))^2)/length(x)
}
# Erro quadrático médio.
cbind(n0 = apply(b0, 1, FUN = eqm),
n1 = apply(b1, 1, FUN = eqm))
#-----------------------------------------------------------------------
# Curvas ajustadas de onde é possível determinar uma banda de confiança.
formula(n0)
formula(n1)
pairs(t(b0))
pairs(t(b1))
par(mfrow = c(1, 2))
plot(Gain ~ A, data = turk0)
apply(b0, MARGIN = 2,
FUN = function(b) {
curve(b[1] + b[2] * x/(b[3] + x),
add = TRUE, col = rgb(0, 0, 1, 0.25), n = 51)
invisible()
})
plot(Gain ~ A, data = turk0)
apply(b1, MARGIN = 2,
FUN = function(b) {
curve(b[1] + b[2] * (1 - 2^(-x/b[3])),
add = TRUE, col = rgb(0, 0, 1, 0.25), n = 51)
invisible()
})
#-----------------------------------------------------------------------
# Estimativa do valor da função em x = 0.2.
p0 <- apply(b0, MARGIN = 2,
FUN = function(b){
b[1] + b[2] * 0.2/(b[3] + 0.2)
})
p1 <- apply(b1, MARGIN = 2,
FUN = function(b){
b[1] + b[2] * (1 - 2^(-0.2/b[3]))
})
# Desvio padrão de \hat{y} em x = 0.2.
cbind(sd(p0), sd(p1))
# Estimativa de "vida" 90%.
q <- 0.9
m0 <- b0[3, ] * q/(1 - q)
m1 <- -log2(1 - q) * b1[3, ]
cbind(n0 = summary(m0), n1 = summary(m1))
# Distribuição do x_{90%}.
par(mfrow = c(1, 2))
plot(density(m0))
rug(m0)
abline(v = coef(n0)["Mei"] * q/(1 - q))
plot(density(m1))
rug(m1)
abline(v = -log2(1 - q) * coef(n1)["Mei"])
layout(1)
plot(Gain ~ A, turk0, xlim = c(0, 2), ylim = c(590, 850))
with(as.list(coef(n0)), {
curve(Int + Ass * A/(Mei + A),
xname = "A", add = TRUE, col = 2)
y90 <- Int + q * Ass
x90 <- Mei * q/(1 - q)
abline(h = y90, v = x90, col = 2, lty = 2)
points(x90, y90, pch = 19, col = 2)
rug(m0, col = 2)
})
with(as.list(coef(n1)), {
curve(Int + Ass * (1 - 2^(-A/Mei)),
xname = "A", add = TRUE, col = 4)
y90 <- Int + q * Ass
x90 <- -log2(1 - q) * Mei
abline(h = y90, v = x90, col = 4, lty = 2)
points(x90, y90, pch = 19, col = 4)
rug(m1, col = 4)
})
#-----------------------------------------------------------------------
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment