ridiculas - Wiki do LEG

Não foi possível enviar o arquivo. Será algum problema com as permissões?

Você está aqui: start » ridiculas

Diferenças

Aqui você vê as diferenças entre duas revisões dessa página.

--- ridiculas [2011/05/18 22:04]
walmes
+++ ridiculas [2011/08/13 07:56]
jcfaria [.Rprofile no linux]
@@ Linha 8: / Linha 8: @@
 ----
-==== Desdobramento de interação em experimento fatorial ====
+==== .Rprofile no linux ====
+<code R>
+#------------------------------------------------------------------------------------------
+#                                                                               por JCFaria
+#------------------------------------------------------------------------------------------
+</code>
+Esse post tem a finalidade de compartilhar algumas coisas que considero importantes na inicialização do R no Linux!
+Muitas das opções importantes, do ponto de vista funcional (não relativos à aparência), podem ser feitas no arquivo .Rprofile.
+Esse arquivo deve ficar localizado no home do usuário (~/.Rprofile) e é um dos primeiros a ser lido quando uma sessão do
+R é iniciada.
+Tenho uma função (bem simples) que uso bastante em meu dia a dia: "cv" para calcular o coef. de variação de uma ANOVA:
+<code R>
+cv <- function(av)
+{
+  if(is.null(av) || !inherits(av, 'aov'))
+    stop('Please, check the parameter!')
+  qmee <- with(av, sum(residuals^2) / df.residual)
+  cv   <- 100 * sqrt(qmee) / mean(av$fitted.values)
+  return(round(cv, 2))
+}
+</code>
+Pois bem, a medida que vamos aumentando nossa intimidade com o R (inevitavelmente) iremos desenvolvendo nossas próprias
+funções (o R foi projetado para isso!).
+<fc #000080>Ai vem o problema: ter que sempre carregar a função quando for usar, o que pode se tornar uma chatisse!</fc>
+Pior ainda, ao limpar o workspace do usuário (.GlobalEnv) elas são removidas e precisam ser recarregadas.
+**Tem como contornar? Sim! De várias formas:**
+**//1. Opção muito pouco prática://**
+<code R>
+oldp <- getwd()
+setwd('/home/jcfaria/dados/r/funcoes/')
+source('cv.r')
+setwd(oldp)
+</code>
+A função "cv" ficará disponível no seu workspace mas será removida com a instrução:
+<code R>
+> rm(list=ls())
+</code>
+muito usada por várias GUIs.
+**//2. Opção "mais" prática://**
+A mesma que a anterior, contudo, a função não deverá se chamar "cv", mas sim ".cv".
+Nesse caso ela permanecerá como um objeto oculto no seu workspace e não será removida com a intrução:
+<code R>
+> rm(list=ls())
+</code>
+Contudo poderá ser removida com a intrução:
+<code R>
+> rm(list=ls(all=TRUE))
+</code>
+**//3. Colocando suas funções em algum ambiente (environment) do R (optei pelo base)://**
+<code R>
+oldp <- getwd()
+setwd('/home/jcfaria/dados/r/funcoes/')
+source('cv.r', local=baseenv())
+setwd(oldp)
+</code>
+Ela não ficará no seu workspace, mas sim no base.
+Como tal, poderá ser usada com qualquer outra função desse pacote.
+**//4. Criando seu próprio ambiente (acho a solução mais elegante)://**
+<code R>
+oldp <- getwd()
+setwd('/home/jcfaria/dados/r/funcoes/')
+.jcf <- new.env()
+source('cv.r', local=.jcf)
+setwd(oldp)
+</code>
+Nesse último caso (**//4//**):
+- O objeto ".jcf" ficará oculto no meu workspace evitando ser deletado com rm(list=ls())
+- Parar acessar a função "cv" será necessário
+<code R>
+> .env$cv
+# ou
+> with(.jcf, cv)
+</code>
+Por exemplo:
+<code R>
+> av <- aov(Sepal.Length ~ Species, data=iris)
+> .jcf$cv(av)
+[1] 8.81
+# ou
+> with(.jcf, cv(av))
+[1] 8.81
+</code>
+Esta forma de carregar funções de forma permanente no R pode ser usado para qualquer outro objeto.
+No Windows bastava usar no /etc/Rprofile.site:
+<code R>
+source('cv', local=TRUE)
+</code>
+que ela ficava disponível no pacote base. Não testei na versão em desenvolvimento (instável) que uso no linux,
+mas deve ainda funcionar.
+----
+==== Como fazer a justaposição de vários data.frames ====
 Temporariamente sem descrição.\\
-palavras-chave: .
+palavras-chave: #merge, #Reduce.
+<code R>
+#------------------------------------------------------------------------------------------
+#                                                                                por Walmes
+#------------------------------------------------------------------------------------------
+id <- 1:30  # número único que identifica os registros
+n <- 20     # número de registros por data.frame
+a1 <- data.frame(id=sample(id, n), v1=rnorm(n))     # resposta 1
+a2 <- data.frame(id=sample(id, n), v2=rpois(n,10))  # resposta 2
+a3 <- data.frame(id=sample(id, n), v3=runif(n))     # resposta 3
+merge(a1, a2, by="id") # justapõe 2 data.frames de cada vez
+a0 <- list(a1, a2, a3) # cria uma lista com todos os data.frames
+Reduce(function(x, y) merge(x, y, by="id"), a0, accumulate=FALSE) # justapõe todos
+#------------------------------------------------------------------------------------------
+# font: http://rwiki.sciviews.org/doku.php?id=tips:data-frames:merge
+#       https://stat.ethz.ch/pipermail/r-help/2008-April/160836.html
+#       http://econometricsense.blogspot.com/2011/01/merging-multiple-data-frames-in-r.html
+#       http://www.youtube.com/watch?v=E4uR5I1uLFM
+#------------------------------------------------------------------------------------------
+</code>
+----
+==== Gráfico de valores observados e curva de valores preditos para modelo linear generalizado ====
+Temporariamente sem descrição.\\
+palavras-chave: #glm, #poisson, #predict.
+<code R>
+#------------------------------------------------------------------------------------------
+#                                                                                por Walmes
+#------------------------------------------------------------------------------------------
+da <- expand.grid(trat=gl(2,1), tempo=1:20)
+y <- rpois(nrow(da), lambda=da$tempo/5)
+g0 <- glm(y~trat/(tempo+I(tempo^2)), data=da, family=poisson)
+new <- expand.grid(trat=gl(2,1), tempo=seq(1,20,l=50))
+new$p0 <- predict(g0, newdata=new, type="response")
+plot(y~tempo, da, col=da$trat)
+with(subset(new, trat=="1"), lines(p0~tempo, col=1))
+with(subset(new, trat=="2"), lines(p0~tempo, col=2))
+#------------------------------------------------------------------------------------------
+</code>
+----
+==== Gráficos de barras com intervalos de confiança para as médias ====
+Temporariamente sem descrição.\\
+palavras-chave: #intervalo, #erro_padrão, #barras.
+<code R>
+#------------------------------------------------------------------------------------------
+#                                                                                por Walmes
+#------------------------------------------------------------------------------------------
+# cmr para colocar intervalos de confiança para as medias num gráfico de barras
+da <- data.frame(trat=gl(5,8))
+da$y <- as.numeric(da$trat)+rnorm(nrow(da))
+m0 <- lm(y~trat, da)
+new <- data.frame(trat=levels(da$trat))
+new$pred <- predict(m0, newdata=new, interval="confidence")
+str(new)
+ylim <- c(0, max(new$pred)*1.05)
+bp <- barplot(new$pred[,1], ylim=ylim)
+arrows(bp, new$pred[,2], bp, new$pred[,3], code=3, angle=90)
+box()
+#------------------------------------------------------------------------------------------
+# outras referências
+browseURL("http://addictedtor.free.fr/graphiques/graphcode.php?graph=54")
+browseURL("http://addictedtor.free.fr/graphiques/RGraphGallery.php?graph=72")
+#------------------------------------------------------------------------------------------
+</code>
+----
+==== Gráficos do R exportados pelo dispositivo tikz ====
+Temporariamente sem descrição.\\
+palavras-chave: #tikz, #sweave, #latex.
+<code>
+\documentclass{article}
+\usepackage{Sweave}
+\usepackage{tikz}
+\SweaveOpts{keep.source=true}
+\title{Usando ti\textit{k}z no Sweave}
+\author{Walmes Zeviani\\ LEG/UFPR}
+\begin{document}
+\maketitle
+C\'{o}digo m\'{i}nimo reproduz\'{i}vel usando \texttt{tikzDevice} para exportar gr\'{a}ficos feitos no R
+para codifica\c{c}\~{a}o ti\textit{k}z.
+{\footnotesize
+<<results=hide>>=
+#------------------------------------------------------------------------------------------
+require(tikzDevice)
+set.seed(2011); x <- rnorm(200)
+tikz("plot.tex", w=5, h=3)
+hist(x, freq=FALSE, ylab="Densidade",
+     main="Histograma de uma amostra de $X \\sim N(\\mu=0, \\sigma^2=1)$")
+curve(dnorm(x), col=2, add=TRUE, lwd=2); rug(x)
+legend("topleft", col=2, lty=1, lwd=2, bty="n",
+       legend="$\\displaystyle \\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\cdot e^{-\\frac{(x-\\mu)^2}{2\\sigma^2}}$")
+box(); dev.off()
+#------------------------------------------------------------------------------------------
+@
+}
+\input{plot.tex}
+\end{document}
+</code>
+----
+==== Gráfico com dois eixos coordenados ====
+Temporariamente sem descrição.\\
+palavras-chave: #série, #eixo.
+<code R>
+#------------------------------------------------------------------------------------------
+#                                                               por Ivan, Benilton e Walmes
+#------------------------------------------------------------------------------------------
+# dados de séries de dados indexadas no tempo (meses)
+lines <- 'meses  temp umidade    rad   chuva
+    Jan    26.49   86.58    795.88  0.36
+    Fev    26.65   88.49    710.24  0.34
+    Mar    27.19   86.16    772.99  0.21
+    Abr    26.28   89.75    574.88  0.67
+    Mai    26.62   89.22    614.02  0.31
+    Jun    26.13   87.83    680.08  0.26
+    Jul    25.83   86.57    675.97  0.15
+    Ago    27.05   83.14    756.44  0.07
+    Set    27.60   83.02    925.57  0.14
+    Out    27.44   85.16    927.71  0.17
+    Nov    26.56   88.18    788.87  0.19
+    Dez    25.87   90.63    703.94  0.33'
+da <- read.table(textConnection(lines), header=TRUE)
+str(da)
+#------------------------------------------------------------------------------------------
+# gráfico de duas séries no mesmo gráfico com um eixo coordenado para cada
+par(mar=c(5,5,4,5))
+plot(da$temp, type="b", pch=15, lwd=1.5,
+     xlab="Meses (Anos 2010)", main="Temperatura e umidade em 2010",
+     ylab=expression(Temperatura~group("(", degree*C, ")")),
+     ylim=c(24,30), xlim=c(1,12), axes=FALSE)
+axis(1, at=1:12, labels=da$meses)
+axis(2)
+par(new=TRUE)
+plot(da$umidade, type="b", pch=14, axes=FALSE, frame=TRUE, ann=FALSE)
+axis(4)
+mtext(text="Umidade (%)", 4, line=3)
+legend("topleft", bty="n", seg.len=3, lty=1, lwd=c(1.5,1), pch=c(15,14),
+       legend=c("Temperatura","Umidade"), merge=TRUE, trace=FALSE)
+#------------------------------------------------------------------------------------------
+# funções sugeridas
+apropos("month")
+demo(plotmath)
+#------------------------------------------------------------------------------------------
+</code>
+----
+==== Análise de dados de proporção usando modelo linear generalizado ====
+Temporariamente sem descrição.\\
+palavras-chave: #binomial, #sucessos, #deviance, #wireframe, #fatorial, #superfície.
+<code R>
+#------------------------------------------------------------------------------------------
+#                                                                                por Walmes
+#------------------------------------------------------------------------------------------
+# dados de número de sementes viáveis de soja
+rend <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/rendimento.txt", header=TRUE)
+rend <- transform(rend, k=factor(K), a=factor(A), bloc=factor(bloc))
+str(rend)
+#------------------------------------------------------------------------------------------
+# ajuste modelo de caselas aos dados assumindo distribuição binomial (link=logit)
+g0 <- glm(cbind(nv, nvi)~bloc+k*a, data=rend, family=binomial)
+#------------------------------------------------------------------------------------------
+# análise de resíduos usual para verificar anomalias
+par(mfrow=c(2,2))
+plot(g0)
+layout(1)
+#------------------------------------------------------------------------------------------
+# quadro de estimativas e quadro de análise de deviance, faz a vez da anova
+summary(g0)
+anova(g0, test="Chisq")
+#------------------------------------------------------------------------------------------
+# obter modelo mais parcimonioso, usar fatores na forma contínua
+g1 <- glm(cbind(nv, nvi)~bloc+K+A+I(K^2)+I(A^2)+K:A, data=rend, family=binomial)
+summary(g1) # comparar deviance residual com grau de liberdade residual
+g1 <- update(g1, formula=.~.-K:A, family=quasibinomial)
+summary(g1)
+anova(g1, test="F")
+#------------------------------------------------------------------------------------------
+# faz a predição dos valores (usando o apenas o bloco 1)
+pred <- with(rend,
+             expand.grid(A=seq(min(A),max(A),l=20),
+                         K=seq(min(K),max(K),l=20),
+                         bloc="1"))
+pred$prob <- predict(g1, newdata=pred, type="response")
+#------------------------------------------------------------------------------------------
+# gráfico
+require(lattice)
+wireframe(prob~A+K, data=pred,
+          zlab=list("Probabilidade de germinação", rot=90),
+          scales=list(arrows=FALSE),
+          screen=list(z=-50, x=-60), drape=TRUE)
+#------------------------------------------------------------------------------------------
+</code>
+----
+==== Desdobramento de interação em experimento fatorial ====
+Temporariamente sem descrição.\\
+palavras-chave: #desdobramento, #comparações, #teste_de_médias, #fatorial.
 <code R>
@@ Linha 81: / Linha 472: @@
 #------------------------------------------------------------------------------------------
-# usando o pacote do ExpDes (versão inglês)
+# usando o ExpDes (inglês) (https://sites.google.com/site/ericbferreira/unifal/downloads-1)
 require(ExpDes)
@@ Linha 157: / Linha 548: @@
               })
-tTukey <- ldply(li, NULL)
+tTukey <- ldply(aux, NULL)
 tTukey$M <- gsub(" ", "", tTukey$M, fixed=TRUE)
 tTukey$trt <- as.factor(as.numeric(as.character(tTukey$trt)))
@@ Linha 163: / Linha 554: @@
 #------------------------------------------------------------------------------------------
-# a cereja em cima da maça como diz o PJ
+# a cereja em cima do bolo como diz o PJ
 barchart(means~trt|.id, data=tTukey, horiz=FALSE, layout=c(3,1),
@@ Linha 226: / Linha 617: @@
 rug(p1)
 abline(v=7, lty=3)
-lines(density(p1), col=1, lwd=2)
+lines(density(p1, from=0, to=10), col=1, lwd=2)
 box()
 #------------------------------------------------------------------------------------------
-</code>
-----
-==== Desdobramento de interação de fatorial duplo em testes de médias ====
-Temporariamente sem descrição.\\
-palavras-chave: #desdobramento, #comparações, #teste_de_médias, #fatorial.
-<code R>
-#------------------------------------------------------------------------------------------
-#                                                                                por Walmes
-#------------------------------------------------------------------------------------------
-# dados
-rend <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/rendimento.txt", header=TRUE)
-rend <- transform(rend, K=factor(K), A=factor(A), bloc=factor(bloc))
-str(rend)
-#------------------------------------------------------------------------------------------
-# ajuste do modelo
-m0 <- aov(rg~bloc+A*K, data=rend)
-summary(m0)
-#------------------------------------------------------------------------------------------
-# checagem
-par(mfrow=c(2,2)); plot(m0); layout(1)
-#------------------------------------------------------------------------------------------
-# desdobrando somas de quadrados para a variação de K dentro de A
-m1 <- aov(rg~bloc+A/K, data=rend)
-summary(m1, split=list("A:K"=list(
-                         "A-37.5"=c(1,4,7,10),
-                         "A-50.0"=c(2,5,8,11),
-                         "A-62.5"=c(3,6,9,12))))
-#------------------------------------------------------------------------------------------
-# para facilitar encontrar as posições pode-se fazer a busca por expessões regulares
-names(coef(m1))[8:19]
-grep("A37.5", names(coef(m1))[8:19])
-grep("A50", names(coef(m1))[8:19])
-grep("A62.5", names(coef(m1))[8:19])
-#------------------------------------------------------------------------------------------
-# usando as expressões regulares vamos desdobrar A dentro de K
-m2 <- aov(rg~bloc+K/A, data=rend)
-summary(m2)
-names(coef(m2))
-#------------------------------------------------------------------------------------------
-# buscando pela expressão regular e criando a lista do desdobramento
-desd <- sapply(paste("K", levels(rend$K), sep=""), simplify=FALSE,
-               function(k){
-                 grep(k, names(coef(m2))[10:19])
-               })
-#------------------------------------------------------------------------------------------
-# decomposição das somas de quadrados
-summary(m2, split=list("K:A"=desd))
-#------------------------------------------------------------------------------------------
-# desdobrando a interação em testes de médias para níveis de K fixando os níveis de A
-require(agricolae)
-with(subset(rend, A=="37.5"),
-     HSD.test(rg, K, DFerror=df.residual(m0), MSerror=deviance(m0)/df.residual(m0)))
-with(subset(rend, A=="50"),
-     HSD.test(rg, K, DFerror=df.residual(m0), MSerror=deviance(m0)/df.residual(m0)))
-with(subset(rend, A=="62.5"),
-     HSD.test(rg, K, DFerror=df.residual(m0), MSerror=deviance(m0)/df.residual(m0)))
-#------------------------------------------------------------------------------------------
-# usando funções para fazer o desdobramento (sapply)
-levels(rend$A)
-sapply(levels(rend$A), simplify=FALSE,
-       function(a){
-         with(subset(rend, A==a),
-              HSD.test(rg, K,
-                       DFerror=df.residual(m0),
-                       MSerror=deviance(m0)/df.residual(m0)))
-       })
-levels(rend$K)
-sapply(levels(rend$K), simplify=FALSE,
-       function(k){
-         with(subset(rend, K==k),
-              HSD.test(rg, A,
-                       DFerror=df.residual(m0),
-                       MSerror=deviance(m0)/df.residual(m0)))
-       })
-#------------------------------------------------------------------------------------------
-# usando o pacote do Eric (https://sites.google.com/site/ericbferreira/unifal/downloads-1)
-require(ExpDes)
-help(package="ExpDes")
-help(fat2.rbd, help_type="html")
-#------------------------------------------------------------------------------------------
-# aplicando a função do Eric fat2.rbd()
-with(rend, fat2.rbd(A, K, bloc, rg, mcomp="tukey", quali=c(TRUE, TRUE)))
-#------------------------------------------------------------------------------------------
 </code>
@@ Linha 1243: / Linha 1522: @@
 hx <- seq(media+2*stder, media+3*stder, .01)
 hy <- dnorm(hx, media, stder)
-n <- length(hy)
+n <-
+ length(hy)
 polygon(c(hx, rev(hx)), c(hy, rep(0, n)), col=2)
 </code>

Diferenças

Navegação

Busca

Ferramentas

QR Code