##======================================================================
## Auxilar function for GLM models por Cesar Taconeli
##======================================================================
envelope <- function(modelo) {
dados <- na.omit(modelo$data)
nsim <- 100
n <- modelo$df.null + 1
r1 <- sort(rstandard(modelo, type = "deviance"))
m1 <- matrix(0, nrow = n, ncol = nsim)
a2 <- simulate(modelo, nsim = nsim)
for (i in 1:nsim) {
dados$y <- a2[, i]
aj <- update(modelo, y ~ ., data = dados)
m1[, i] <- sort(rstandard(aj, type = "deviance"))
}
li <- apply(m1, 1, quantile, 0.025)
m <- apply(m1, 1, quantile, 0.5)
ls <- apply(m1, 1, quantile, 0.975)
quantis <- qnorm((1:n - 0.5)/n)
plot(rep(quantis, 2), c(li, ls), type = "n",
xlab = "Percentis teóricos",
ylab = "Resíduos")
## title('Gráfico Normal de Probabilidades')
lines(quantis, li, type = "l")
lines(quantis, m, type = "l", lty = 2)
lines(quantis, ls, type = "l")
points(quantis, r1, pch = 16, cex = 0.75)
}
CHosmer <- function(modelo, g) {
respostas <- modelo$y
## Dataframe com probabilidades estimadas e respostas para cada
## indivíduo.
preditos <- predict(modelo, type = "response")
## Ordenando as linhas do dataframe da menor para a maior
## probabilidade estimada.
dpred <- data.frame(preditos, respostas)
## Calculando os quantis para as probabilidades estimadas, para
## posterior formação dos grupos.
dpred <- dpred[order(dpred[, 1]), ]
## Formando g grupos, de tamanhos (aproximadamente) iguais, com
## probabilidades estimadas semelhantes.
cortes <- quantile(dpred[, 1], probs = seq(0, 1, 1/g),
include.lowest = TRUE)
c1 <- cut(dpred[, 1], breaks = cortes, include.lowest = T)
## Obs é um vetor com o número observado de respostas em cada um dos
## g grupos.
Obs <- tapply(dpred[, 2], c1, sum)
## pi é um vetor com as médias das probabilidades estimadas em cada
## um dos g grupos.
pi <- tapply(dpred[, 1], c1, mean)
## n é um vetor com os tamanhos de amostras em cada grupo.
n <- tapply(dpred[, 1], c1, length)
## Estatística do teste de qualidade proposto por Hosmer e Lemeshow e
Cchap <- sum(((Obs - n * pi)^2)/(n * pi * (1 - pi)))
## p-valor correspondente.
pv <- 1 - pchisq(Cchap, g - 2)
return(list(Cchap = Cchap, pvalue = pv))
}
diag.binom <- function(fit.model){
X <- model.matrix(fit.model)
n <- nrow(X)
p <- ncol(X)
w <- fit.model$weights
W <- diag(w)
H <- solve(t(X)%*%W%*%X)
H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
h <- diag(H)
ts <- resid(fit.model,type="pearson")/sqrt(1-h)
td <- resid(fit.model,type="deviance")/sqrt(1-h)
di <- (h/(1-h))*(ts^2)
a <- max(td)
b <- min(td)
par(mfrow=c(1,3))
plot(td,xlab="Índice", ylab="Resíduo Componente do Desvio",
ylim=c(b-1,a+1), pch=20)
title(sub="(a)")
abline(h=c(2, -2),lty=2)
plot(fitted(fit.model), td,xlab="Valor Ajustado",
ylab="Resíduo Componente do Desvio", pch=20)
title(sub="(b)")
abline(h=c(-2, 0, 2), lty=2)
envelope(fit.model)
title(sub="(c)")
}
diag.pois <- function(fit.model, da, ni){
X <- model.matrix(fit.model)
n <- nrow(X)
p <- ncol(X)
w <- fit.model$weights
W <- diag(w)
H <- solve(t(X)%*%W%*%X)
H <- sqrt(W)%*%X%*%H%*%t(X)%*%sqrt(W)
h <- diag(H)
ts <- resid(fit.model,type="pearson")/sqrt(1-h)
td <- resid(fit.model,type="deviance")/sqrt(1-h)
di <- (h/(1-h))*(ts^2)
par(mfrow=c(3,2))
a <- min(td)
b <- max(td)
plot(fitted(fit.model), h,xlab="Valor Ajustado", ylab="Medida h",
pch=16)
## identify(fitted(fit.model), h, n=1)
##
envelope(m3o)
## identify(di, n=1)
##
plot(td,xlab="Índice", ylab="Resíduo Componente do Desvio",
ylim=c(a-1,b+1), pch=16)
abline(2,0,lty=2)
abline(-2,0,lty=2)
## identify(td, n=1)
##
w <- fit.model$weights
eta <- predict(fit.model)
z <- eta + resid(fit.model, type="pearson")/sqrt(w)
plot(predict(fit.model),z,xlab="Preditor Linear",
ylab="Variavel z", pch=16)
lines(smooth.spline(predict(fit.model), z, df=2))
hatv <- ((influence.measures(m3o)$infmat[,15]))
dicook <- ((influence.measures(m3o)$infmat[,14]))
plot(hatv, ylab="Alavancagem", type="n")
grid()
lines(x=1:nrow(da), y=hatv, type="h", ylab="Alavancagem")
points(x=1:nrow(da), y=hatv, pch=20)
text(x = 45 - ni, hatv[45], "45", cex=0.95)
##
plot(dicook, ylab="Dist. de Cook", type="n")
grid()
lines(x=1:nrow(da), y=dicook, type="h")
points(x=1:nrow(da), y=dicook, pch=20)
text(x = 69 - ni, dicook[69], "69", cex=0.95)
text(x = 45 - ni, dicook[45], "45", cex=0.95)
par(mfrow=c(1,1))
}
##======================================================================
## Lattice functions por Walmes Zeviane
## Funções desenvolvidas para matérias do blog 'Ridículas <-
## dicas curtas sobre R' www.ridiculas.wordpress.com, por
## Walmes Zeviani envie sugestões para walmes<at>ufpr.br
##======================================================================
library(lattice)
library(latticeExtra)
##----------------------------------------------------------------------
## painel para gráficos com contornos
panel.3d.contour <- function(x, y, z, rot.mat, distance, type = "on",
nlevels = 20, zlim.scaled, col.contour = 1, ...) {
clines <- contourLines(x, y, matrix(z, nrow = length(x),
byrow = TRUE), nlevels = nlevels)
if (any(type %in% c("bottom"))) {
for (ll in clines) {
n <- ltransform3dto3d(rbind(ll$x, ll$y, zlim.scaled[1]),
rot.mat, distance)
panel.lines(n[1, ], n[2, ], col = col.contour, lty = 1,
lwd = 1)
}
}
panel.3dwire(x, y, z, rot.mat, distance, zlim.scaled = zlim.scaled,
...)
if (any(type %in% c("on"))) {
for (ll in clines) {
n <- ltransform3dto3d(rbind(ll$x, ll$y, ll$level),
rot.mat, distance)
panel.lines(n[1, ], n[2, ], col = col.contour, lty = 1,
lwd = 1)
}
}
if (any(type %in% c("top"))) {
for (ll in clines) {
n <- ltransform3dto3d(rbind(ll$x, ll$y, zlim.scaled[2]),
rot.mat, distance)
panel.lines(n[1, ], n[2, ], col = col.contour, lty = 1,
lwd = 1)
}
}
}
## Exemplo de uso
da <- expand.grid(x = seq(0, 1, l = 20), y = seq(0, 1, l = 20))
da$z <- with(da, 10 + x + 0.5 * y + 2 * x * y - 1 * x^2 - 1 * y^2)
wireframe(z ~ x + y, da, zlim = c(8, 12), drape = TRUE, col = "white",
col.contour = "red", panel.3d.wireframe = "panel.3d.contour",
type = c("on", "bottom"), screen = list(z = 40, x = -75))
##----------------------------------------------------------------------
## painéis para fazer curvas de regressão com bandas de
## confiança
prepanel.ciH <- function(x, y, ly, uy, subscripts, ...) {
x <- as.numeric(x)
ly <- as.numeric(ly[subscripts])
uy <- as.numeric(uy[subscripts])
list(ylim = range(uy, ly, finite = TRUE), xlim = range(x))
}
panel.ciH <- function(x, y, ly, uy, subscripts, ...) {
y <- as.numeric(y)
x <- as.numeric(x)
or <- order(x)
ly <- as.numeric(ly[subscripts])
uy <- as.numeric(uy[subscripts])
panel.polygon(c(x[or], x[rev(or)]), c(ly[or], uy[rev(or)]),
col = 1, alpha = 0.15, border = NA)
panel.lines(x[or], ly[or], lty = 3, lwd = 0.5, col = 1)
panel.lines(x[or], uy[or], lty = 3, lwd = 0.5, col = 1)
panel.xyplot(x, y, subscripts = subscripts, ...)
}
## ## Exemplo de uso
## da <- expand.grid(A = gl(2, 4), B = 0:10)
## da$y <- with(da, rnorm(A, mean = as.numeric(A) * B, sd = 2))
## m0 <- lm(y ~ A * B, da)
## da <- cbind(da, predict(m0, interval = "confidence"))
## str(da)
##
## ## Sem grupos
## xyplot(y ~ B | A, da) +
## as.layer(
## xyplot(fit ~ B | A,
## da,
## type = "l",
## ly = da$lwr,
## uy = da$upr,
## prepanel = prepanel.ciH,
## panel = panel.ciH))
##
## ## Com grupos
## xyplot(y ~ B, groups = A, da) +
## as.layer(
## xyplot(fit ~ B,
## groups = A,
## da,
## type = "l",
## ly = da$lwr,
## uy = da$upr,
## prepanel = prepanel.ciH,
## panel = panel.superpose,
## panel.groups = panel.ciH))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.