Archive

Posts Tagged ‘substitute’

Animações em documentos PDF produzidos com Sweave

Animação que representa as probabilidades obtidas pela tabela da distribuição normal padrão.

O pacote do animate do \LaTeX, desenvolvido por Alexander Grahn, permite que sejam inseridas animações em documentos PDF. Esses arquivos devem ser visualizados em leitores de PDF da Adobe. No Linux podemos usar o acroread.

O procedimento é simples. Vou fazer isso em um documento Sweave. Primeiro confeccionamos diversas figuras que formam uma seqüência de ações. Tais figuras podem ser feitas com controladores de fluxo como a função for() e similares. Inclusive, o pacote animation do R, desenvolvido por Yihui Xie, oferece implementações para confecção de animações. As figuras são salvas em um diretório. Para adicioná-las ao documento na forma de animação (com botões de play, stop) basta usar o comando \animategraphics{}.

O arquivo Sweave abaixo (*.Rnw) gera um arquivo PDF com a tabela da distribuição normal padrão e uma animação que representa a área integrada correspondente aos valores da tabela. Veja o PDF. Para gerar gifs animados como R veja o artigo do Mark Heckmann no R you ready?. Até a próxima ridícula.

\documentclass[a4paper]{article}
\usepackage[brazil]{babel}
\usepackage{Sweave}
\usepackage{animate}
\thispagestyle{empty}

\begin{document}

\begin{center}
{\Large Tabela da distribui\c{c}\~{a}o normal padr\~{a}o}\\ \vspace{2ex}
{\large Walmes Zeviani -- LEG/UFPR}  \\
\end{center}

<<echo=false, results=hide>>=
#-----------------------------------------------------------------------------
dir.create("exemplos")
png(file="exemplos/qnorm%1d.png", width=500, height=250)
par(mar=c(4,4,1,1))
for(q in seq(0, 4,l=100)){
  curve(dnorm(x, 0, 1), -5, 5, ylab="f(z)", xlab="z")
  x <- seq(0, q, by=0.01)
  fx <- dnorm(x, 0, 1)
  polygon(c(x, rev(x)),
          c(fx, rep(0, length(fx))),
          col="gray90")
  abline(v=0, lty=2)
  Pr <- round(pnorm(q, 0, 1)-0.5, digits=3)
  qq <- round(q, digits=3)
  legend("topleft", bty="n", fill="gray90",
         legend=substitute(P(0<~Z<=~q)==Pr, list(q=qq, Pr=Pr)))
}
dev.off()
#-----------------------------------------------------------------------------
@ 

<<echo=false, results=hide>>=
#-----------------------------------------------------------------------------
require(xtable)
options(OutDec=",")
q <- seq(0,3.99,by=0.01)
p <- pnorm(q)-0.5
m <- matrix(p, byrow=TRUE, ncol=10)
rownames(m) <- gsub("\\.", ",", formatC(seq(0,3.9,0.1), dig=1, format="f"))
colnames(m) <- 0:9/100
#-----------------------------------------------------------------------------
@

\begin{center}
\animategraphics[controls, loop, width=0.75\textwidth]{10}{exemplos/qnorm}{1}{100}
\end{center}

\begin{center}
\small\addtolength{\tabcolsep}{-3pt}
{\footnotesize
<<echo=false, results=tex>>=
#-----------------------------------------------------------------------------
print(xtable(m, digits=5), floating=FALSE)
#-----------------------------------------------------------------------------
@ 
}
\end{center}
\end{document}

Expressões, bandas de confiança e legenda

Nesse post vou apresentar como confeccionar um gráfico para representar o ajuste de um modelo de regressão. Esse gráfico vai conter elementos básicos de um ajuste, ou seja, valores observados, valores preditos, bandas de confiança, legenda e equação estimada (veja figura abaixo), que são elementos que normalmente aparecem em gráficos de artigos científicos e gráficos gerados por outros aplicativos estatísticos (menos interessantes que o R, lógico!).

Valores observados e preditos pelo modelo de regressão com o intervalo de confiança (95%).

Vamos usar os dados contidos no data.frame “anscombe”. Esse conjunto de dados é usado em livros de regressão para exemplificar a utilidade das medidas de diagnóstico de resíduos. Em outra oportunidade vou postar algo sobre isso.

No script abaixo, na linha [10] foi ajustado o modelo de equação da reta

y_{1i} = \beta_0 + \beta_1 \cdot x_{1i} + \epsilon_i.

aos dados. O resumo desse ajuste é apresentado em [11]. Foi feito a predição da resposta segundo a equação acima e as estimativas dos parâmetros no passo [16:17]. Usamos um intervalo de 50 pontos para obtermos linhas suaves. As estimativas dos parâmetros e R² foram arredondadas e armazenadas em um vetor [19].

No último bloco do script, trocamos a opção para que os decimais sejam representados por vírgula tanto no console quanto nos gráficos [24], pois na língua portuguesa o separador decimal é a virgula. Em [25] fazemos o diagrama de dispersão dos dados observados, adicionando as retas do ajuste e das bandas de confiança em [26]. Em [27:30] é colocada a legenda ao gráfico no canto inferior direito, identificando com ponto de linhas os valores observados, valores preditos e as bandas de confiança de 95% para o valor predito. Finalmente, em [31:33] é adicionado à equação ajustada ao gráfico. Usamos a função locator() para que o usuário clique no gráfico para adicionar a equação. A função substitute() foi usada para automaticamente substituir os valores dos parâmetros por suas estimativas na equação.

#-----------------------------------------------------------------------------
# dados

data(anscombe)
str(anscombe)

#-----------------------------------------------------------------------------
# ajuste do modelo

m0 <- lm(y1~x1, data=anscombe)
summary(m0)

#-----------------------------------------------------------------------------
# fazendo a predição intervalar num grid regular mais fino

pred <- data.frame(x1=seq(min(anscombe$x1), max(anscombe$x1), length=50))
pred$y1 <- predict(m0, newdata=pred, interval="confidence")
str(pred)
cf <- format(c(coef(m0), summary(m0)$r.squared), digits=4)

#-----------------------------------------------------------------------------
# fazendo o gráfico, descomente para salvar a figura em png

#png("f001.png", w=500, h=300); par(mar=c(5.1,4.1,2.1,2.1))
options(OutDec=",")
plot(y1~x1, data=anscombe, xlab=expression(x[1]), ylab=expression(y[1]))
with(pred, matlines(x1, y1, col=c(1,2,2), lty=c(1,2,2)))
legend("bottomright",
       legend=c("valores observados", "valores preditos",
         "intervalo de confiança (95%)"),
       lty=c(NA,1,2), col=c(1,1,2), pch=c(1,NA,NA), bty="n")
text(locator(n=1), # ao salvar a figura use x=7.5, y=10 no lugar da locator()
     label=substitute(hat(y)[1]==a+b%.%x[1]~~~(R^2==c),
       list(a=cf[1], b=cf[2], c=cf[3]))) # clicar no gráfico
#dev.off()

#-----------------------------------------------------------------------------