Não foi possível enviar o arquivo. Será algum problema com as permissões?
Diferenças
Aqui você vê as diferenças entre duas revisões dessa página.
Ambos lados da revisão anterior Revisão anterior Próxima revisão | Revisão anterior | ||
ridiculas [2011/05/19 18:55] walmes [section 2] |
ridiculas [2012/10/25 15:45] (atual) walmes [section 8] |
||
---|---|---|---|
Linha 1: | Linha 1: | ||
- | ===== Ridículas - dicas curtas sobre R ===== | + | ===== R-idículas - dicas curtas sobre R ===== |
{{:nuvem.png?550 |}} | {{:nuvem.png?550 |}} | ||
- | **//Ridículas//** é a página do LEG dedicada à fornecer //dicas curtas// sobre R, e.g. condução de análises, operação com dados e confecção de gráficos. As dicas estão organizadas pelo título, seguido de descrição, palavras-chave e CMR (código mínimo reproduzível). Se você deseja contribuir com a nossa página de Ridículas, envie e-mail para ''walmes@ufpr.br''. | + | **//R-idículas//** é a página do LEG dedicada à fornecer //dicas curtas// sobre R, e.g. entre outras: otimização do ambiente, condução de análises, operação com dados e confecção de gráficos. As dicas estão organizadas pelo título, seguido de descrição, palavras-chave e CMR (código mínimo reproduzível). Se você deseja contribuir com a nossa página de R-idículas, envie e-mail para ''walmes@ufpr.br''. |
+ | |||
+ | ---- | ||
+ | |||
+ | ==== Controlando gráficamente parâmetros de distribuições de probabilidade ==== | ||
+ | |||
+ | Temporariamente sem descrição.\\ | ||
+ | palavras-chave: #rpanel. | ||
+ | |||
+ | <code R> | ||
+ | #------------------------------------------------------------------------------------------ | ||
+ | # por Walmes | ||
+ | #------------------------------------------------------------------------------------------ | ||
+ | |||
+ | require(rpanel) | ||
+ | |||
+ | norm.panel <- function(panel){ | ||
+ | ##------------------------------------------------------------------ | ||
+ | ## panel$interval: vetor com domínio de plotagem da função | ||
+ | ## penel$...: serão parâmetros da distribuição de probabilidades | ||
+ | curve(dnorm(x, mean=panel$mean, sd=panel$sd), | ||
+ | from=panel$interval[1], to=panel$interval[2]) | ||
+ | panel | ||
+ | ##------------------------------------------------------------------ | ||
+ | } | ||
+ | |||
+ | # passar os argumentos que serão fixos, abre a janelinha | ||
+ | panel <- rp.control(interval=c(-4,4)) | ||
+ | |||
+ | # controla a média | ||
+ | rp.slider(panel, mean, -4, 4, initval=0, showvalue=TRUE, action=norm.panel) | ||
+ | |||
+ | # controla o desvio-padrão | ||
+ | rp.slider(panel, sd, 0.001, 10, initval=1, showvalue=TRUE, action=norm.panel) | ||
+ | |||
+ | #------------------------------------------------------------------------------------------ | ||
+ | </code> | ||
+ | |||
+ | ---- | ||
+ | |||
+ | |||
+ | ==== .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! | ||
+ | Pior ainda, ao limpar o workspace do usuário (.GlobalEnv) ela é removida e precisa ser recarregada. Imagine o contexto da necessidade de várias funções, que é o que geralmente acontece com o passar do tempo ...</fc> | ||
+ | |||
+ | <fc #000080>**Tem como contornar? Sim! De várias formas:**</fc> | ||
+ | |||
+ | **//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 meu 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 meu 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 meu 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" é necessário: | ||
+ | |||
+ | <code R> | ||
+ | > .jcf$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. | ||
+ | |||
+ | No Windows não testei na versão em desenvolvimento (instável) que uso no Linux, | ||
+ | mas deve funcionar, pois usei por muitos anos em várias versões. | ||
+ | |||
+ | Abaixo meu .Rprofile: | ||
+ | <code R> | ||
+ | ## José Cláudio Faria/UESC/DCET | ||
+ | |||
+ | ##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬ | ||
+ | ## General options | ||
+ | library(utils) | ||
+ | options(list(repos='http://cran-r.c3sl.ufpr.br/', | ||
+ | width=80, | ||
+ | editor='vim', | ||
+ | browser='chromium')) | ||
+ | |||
+ | ##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬ | ||
+ | ## Package debug: options | ||
+ | options(list(debug.height=10, | ||
+ | debug.width=125, | ||
+ | debug.screen.pos='-1680+00' # dois monitores! | ||
+ | )) | ||
+ | |||
+ | ##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬ | ||
+ | ## Funtions: myself | ||
+ | oldp <- getwd() | ||
+ | setwd('/home/jcfaria/dados/r/funcoes/') | ||
+ | #source('fitreg.r', local=baseenv()) | ||
+ | #source('fitregl.r', local=baseenv()) | ||
+ | #source('fitrlcor.r', local=baseenv()) | ||
+ | #source('plotreg.r', local=baseenv()) | ||
+ | #source('plotrl.r', local=baseenv()) | ||
+ | #source('plotlm3d.r', local=baseenv()) | ||
+ | #source('cv.r', local=baseenv()) | ||
+ | #source('rs.r', local=baseenv()) | ||
+ | #source('sourcedir.r', local=baseenv()) | ||
+ | |||
+ | .jcf <- new.env() | ||
+ | source('fitreg.r', local=.jcf) | ||
+ | source('fitregl.r', local=.jcf) | ||
+ | source('fitrlcor.r', local=.jcf) | ||
+ | source('plotreg.r', local=.jcf) | ||
+ | source('plotrl.r', local=.jcf) | ||
+ | source('plotlm3d.r', local=.jcf) | ||
+ | source('cv.r', local=.jcf) | ||
+ | source('rs.r', local=.jcf) | ||
+ | source('sourcedir.r', local=.jcf) | ||
+ | #attach(.jcf) | ||
+ | |||
+ | ##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬ | ||
+ | ## Removing variables | ||
+ | setwd(oldp) | ||
+ | rm(oldp) | ||
+ | |||
+ | ##¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬ | ||
+ | ## Packages: loading | ||
+ | #library(fdth) | ||
+ | #library(ScottKnott) | ||
+ | #library(bpca) | ||
+ | #library(TukeyC) | ||
+ | #library(reshape) | ||
+ | #library(debug) | ||
+ | </code> | ||
+ | |||
+ | ---- | ||
+ | |||
+ | |||
+ | ==== Como fazer a justaposição de vários data.frames ==== | ||
+ | |||
+ | Temporariamente sem descrição.\\ | ||
+ | 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 | ||
+ | |||
+ | #------------------------------------------------------------------------------------------ | ||
+ | |||
+ | mergeAll <- function(..., by="date", all=TRUE){ | ||
+ | dotArgs <- list(...) | ||
+ | Reduce( | ||
+ | function(x, y){ | ||
+ | merge(x, y, by = by, all = all, suffixes=paste(".", names(dotArgs), sep = "")) | ||
+ | }, | ||
+ | dotArgs) | ||
+ | } | ||
+ | |||
+ | mergeAll(a = a1, b = a2, c = a3, by="id") | ||
+ | |||
+ | str(.Last.value) | ||
+ | |||
+ | You also might be able to set it up to capture names without you having to put "a = a" etc. using substitute. | ||
+ | |||
+ | # http://r.789695.n4.nabble.com/merge-multiple-data-frames-td4331089.html | ||
+ | # ver reshape::merge_all() | ||
+ | #------------------------------------------------------------------------------------------ | ||
+ | </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) | ||
+ | |||
+ | #------------------------------------------------------------------------------------------ | ||
+ | # para enviar para outra pessoa pode-se enviar o resultado do comando dput(), assim | ||
+ | |||
+ | dput(da) | ||
+ | |||
+ | # cole na mensagem o resulado que aparece no console ao executar essa linha | ||
+ | # assim a pessoa pode ler fazendo assim | ||
+ | |||
+ | da <- | ||
+ | structure(list(meses = structure(c(5L, 4L, 9L, 1L, 8L, 7L, 6L, | ||
+ | 2L, 12L, 11L, 10L, 3L), .Label = c("Abr", "Ago", "Dez", "Fev", | ||
+ | "Jan", "Jul", "Jun", "Mai", "Mar", "Nov", "Out", "Set"), class = "factor"), | ||
+ | temp = c(26.49, 26.65, 27.19, 26.28, 26.62, 26.13, 25.83, | ||
+ | 27.05, 27.6, 27.44, 26.56, 25.87), umidade = c(86.58, 88.49, | ||
+ | 86.16, 89.75, 89.22, 87.83, 86.57, 83.14, 83.02, 85.16, 88.18, | ||
+ | 90.63), rad = c(795.88, 710.24, 772.99, 574.88, 614.02, 680.08, | ||
+ | 675.97, 756.44, 925.57, 927.71, 788.87, 703.94), chuva = c(0.36, | ||
+ | 0.34, 0.21, 0.67, 0.31, 0.26, 0.15, 0.07, 0.14, 0.17, 0.19, | ||
+ | 0.33)), .Names = c("meses", "temp", "umidade", "rad", "chuva" | ||
+ | ), class = "data.frame", row.names = c(NA, -12L)) | ||
+ | |||
+ | #------------------------------------------------------------------------------------------ | ||
+ | # 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> | ||
---- | ---- | ||
Linha 225: | Linha 756: | ||
rug(p1) | rug(p1) | ||
abline(v=7, lty=3) | abline(v=7, lty=3) | ||
- | lines(density(p1), col=1, lwd=2) | + | lines(density(p1, from=0, to=10), col=1, lwd=2) |
box() | box() | ||
Linha 233: | Linha 764: | ||
---- | ---- | ||
- | ==== Desdobramento de interação usando a multcomp::ghlt() ==== | + | ==== Desdobramento de interação usando a multcomp::glht() ==== |
Temporariamente sem descrição.\\ | Temporariamente sem descrição.\\ | ||
Linha 1130: | Linha 1661: | ||
hx <- seq(media+2*stder, media+3*stder, .01) | hx <- seq(media+2*stder, media+3*stder, .01) | ||
hy <- dnorm(hx, media, stder) | hy <- dnorm(hx, media, stder) | ||
- | n <- length(hy) | + | n <- |
+ | length(hy) | ||
polygon(c(hx, rev(hx)), c(hy, rep(0, n)), col=2) | polygon(c(hx, rev(hx)), c(hy, rep(0, n)), col=2) | ||
</code> | </code> | ||