Outubro 2013 - R-br

Variável Dummy na regressão
by Edson Lira 30 Out '13

30 Out '13

Caros amigos, estou tentando fazer um estudo de regressão com duas variáveis preditoras. Dados e comandos abaixo: mus<-c(2,6,8,3,2,7,9,8,4,6) tr<-c("e","m","e","m","e","e","m","m","e","e") trh<-c(2.9,3.0,4.8,1.8,2.9,4.9,4.2,4.8,4.4,4.5) tipo<-ifelse(tr=="e",1,0) #codificar a variável tr, 1=eletrico 2=mecânico b1<-data.frame(mus,tipo,trh) Quero fazer a análise usando a variável tipo. São dois modelos: um quando o e=1 e o outro e=0. m1<-lm(trh~mus+tipo=0, data=b1) m2<-lm(trh~mus+tipo=1, data=b1) Não estou conseguindo implementar, alguém poderia dar alguma sugestão? [ ]'s. Edson Lira Estatístico Manaus-Amazonas

4 4

Re: [R-br] erro na validação cruzada com pacote geoR
by Éder Comunello 30 Out '13

30 Out '13

Senhores, bom dia! Agradeço novamente ao Prof. Paulo pela atenção dispensada. Hélio, quanto ao script que eu tinha ficado de olhar, vou colocar minhas considerações abaixo. Procurei não mudar muito, mas acabei alterando um pouco. Os comentários estão no corpo do script. Lembrando sempre que também estou aprendendo sobre o assunto, então pode ter muita coisa errada ou mal interpretada aí. Mas acho que pode te ajudar. Outra coisa é que não tenho muito ideia do que os dados são, como foram obtidos e as condições da área de onde vieram, o que pra mim é fundamental pra análise. Pelos indícios nos dados (localização, estrutura), vou arriscar que são dados da produção de café ou citros e que o terreno está numa pendente da direita pra esquerda, com maior produção na parte baixa. O que chamou a atenção é o formato da área que é bem particular e não parece muito 'natural'. Mas como o foco principal é o uso do R/geoR, vou tentar me deter nisso. Use por sua própria conta e risco! :D ### <BEGIN> require(geoR) #;require(tcltk2) setwd("C:/LAB/RGIS/trend"); dir() par.ori <- par(no.readonly = TRUE) url1 <- 'https://dl.dropboxusercontent.com/u/117618178/helio/dados.txt' url2 <- 'https://dl.dropboxusercontent.com/u/117618178/helio/recorte.txt' ### Como criar uma pasta pública no Dropbox... browseURL('https://www.dropbox.com/help/16/pt_BR') ### 1. Carregando dados ------------------------------------------------------ dados <- read.geodata(url1, coords.col=1:2, data.col=3, sep="", header=T) dados$borders <- read.table(url2, header=T) summary(dados); length(dados$data) ls(dados);ls.str(dados) ### 2. Visualizando dados ---------------------------------------------------- ### Visão geral par(mfrow=c(1,1)); par()$mfrow points(dados) points(dados, pt.div="quart") points(dados, pt.div="quint") ### tem pelos três sub-áreas bem marcadas... ### Verificando tendências com o plot.geodata -------------------------------- ### Antes de partir pros variogramas, use o plot.geodata para pré-avaliar tendências ### Bom para esclarecer opções de trend, pois notei certa confusão (cte,1st,2nd) ### Recomendo ler o item 'trend' na ajuda da função plot.geodata ?plot.geodata plot(dados, lowess=T) plot(dados, lowess=T, trend='cte') ### mesmo que anterior ### linhas de comando anteriores são equivalentes: trend='cte' é 'default' da geoR! ### nos gráf. 2 e 3, plota-se as coord vs. dados. Observar as distorções na linha 'lowess' plot(dados, lowess=T, trend=~coords) plot(dados, lowess=T, trend='1st') plot(dados, lowess=T, trend=~dados$coords[,1]+dados$coords[,2]) ### as três linhas anteriores são equivalentes ### o uso da polinomial de primeira ordem melhora o comportamento do histograma ### mas a linha de lowess continua com comportamento anômalo plot(dados, lowess=T, trend='2nd') plot(dados, lowess=T, trend=~I(coords)+I(coords^2)+I(coords[,1]*coords[,2])) ### as duas linhas anteriores são equivalentes (polinomal de segunda ordem)! ### não aparenta agregar nenhuma melhoria ### Mais sobre modelos de tendência com trend.spatial ------------------------ ### Pra melhorar o entendimento dos modelos, dá pra dar uma olhada nas matrizes ?trend.spatial head(unclass(trend.spatial('cte', dados))) ### constante head(unclass(trend.spatial(~coords, dados))) ### polinomal 1a ordem, usando coords head(unclass(trend.spatial('1st', dados))) ### mesmo que anterior head(unclass(trend.spatial('2nd', dados))) ### polinomal 2a ordem, usando coords head(unclass(trend.spatial(~coords[,1], dados))) ### em função de X head(unclass(trend.spatial(~coords[,2], dados))) ### em função de Y head(unclass(trend.spatial(~coords+I(coords[,1]^2), dados))) ### coords + X**2 head(unclass(trend.spatial(~I(coords)+I(coords^2), dados))) head(unclass(trend.spatial(~I(coords)+I(coords^2)+I(coords[,1]*coords[,2]), dados))) ### mesmo que 2nd head(unclass(trend.spatial(~I(coords)+I(coords[,1]^2), dados))) head(unclass(trend.spatial(~I(coords^2), dados))) head(unclass(trend.spatial(~coords[,1]*coords[,2], dados))) head(unclass(trend.spatial(~I(coords[,1]*coords[,2]), dados))) head(unclass(trend.spatial(~coords[,1]+I(coords[,1]^2)+coords[,2], dados))) ### Separando o drift com base na matriz trend.spatial ----------------------- trend.matrix <- unclass(trend.spatial(trend='1st', geodata=dados)); nrow(trend.matrix) trend.lm <- lm(dados$data ~ trend.matrix) trend.res <- lm(dados$data ~ trend.matrix)$residuals trend.fit <- lm(dados$data ~ trend.matrix)$fitted.values ### 'Empacotando' os resíduos em um geodata----------------------------------- dados.res <- dados names(trend.res) <- NULL; trend.res dados.res$data <- trend.res; dados.res$data str(dados.res) ### Pode usar dados.res pra fazer a análise variográfica no lugar de dados, ### caso tenha problemas com krige.conv()!!! ### Embora o uso da polinomial de 1a ordem tenha representado alguma melhoria, ### eu não fiquei muito entusiamado por conta do 'lowess' apresentado ### Talvez fosse o caso de testar covariáveis. ### Uma fácil de obter e que pode ter relação é a altitude... ### Mas vamos investigar o efeito nos variogramas ### 3. Distancia máxima e parâmetros para semivariograma --------------------- hmax <- summary(dados)[[3]][[2]]*.8 ; hmax ### 80% da máxima distância uvec <- 12; pairs <- 20 ### pairs.min ### 4. Variogramas ----------------------------------------------------------- ### Não precisaria rodar todos, afinal como visto, tá repetindo modelos!!! ### v0:média constante; v1:1st; v2:2nd; v3:cte; v4:coords v0 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical") v1 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical", trend="1st") v2 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical", trend="2nd") v3 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical", trend="cte") v4 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical", trend=~coords) ### envelopes v0.env <- variog.mc.env(dados, obj.var = v0) v1.env <- variog.mc.env(dados, obj.var = v1) v2.env <- variog.mc.env(dados, obj.var = v2) v3.env <- variog.mc.env(dados, obj.var = v3) v4.env <- variog.mc.env(dados, obj.var = v4) ### 5. Visualizando variogramas ---------------------------------------------- par(par.ori) par(mfrow=c(2,3)) points(dados, pt.div="quartiles",xlab="Leste", ylab="Norte") plot(v0, xlab="Dist. (m)", ylim=c(0,summary(v0$v)[6]*1.1), env=v0.env, main="Isotrópico com tendência") text(v0$u,v0$v,round(v0$n,1),pos=1);text(v0$u,v0$v,round(v0$u,1),pos=3);lines(v0$u,v0$v);abline(h=mean(v0$v)) plot(v1, xlab="Dist. (m)", ylim=c(0,summary(v1$v)[6]*1.1), env=v1.env, main="Isotrópico trend 1st") text(v1$u,v1$v,round(v1$n,1),pos=1);text(v1$u,v1$v,round(v1$u,1),pos=3);lines(v1$u,v1$v);abline(h=mean(v1$v)) plot(v2, xlab="Dist. (m)", ylim=c(0,summary(v2$v)[6]*1.1), env=v2.env, main="Isotrópico trend 2nd") text(v2$u,v2$v,round(v2$n,1),pos=1);text(v2$u,v2$v,round(v2$u,1),pos=3);lines(v2$u,v2$v);abline(h=mean(v2$v)) plot(v3, xlab="Dist. (m)", ylim=c(0,summary(v3$v)[6]*1.1), env=v3.env, main="Isotrópico trend cte") text(v3$u,v3$v,round(v3$n,1),pos=1);text(v3$u,v3$v,round(v3$u,1),pos=3);lines(v3$u,v3$v);abline(h=mean(v3$v)) plot(v4, xlab="Dist. (m)", ylim=c(0,summary(v4$v)[6]*1.1), env=v4.env, main="Isotrópico trend coords") text(v4$u,v4$v,round(v4$n,1),pos=1);text(v4$u,v4$v,round(v4$u,1),pos=3);lines(v4$u,v4$v);abline(h=mean(v4$v)) ### Analisando os variogramas ------------------------------------------------ ### Realmente o comportamento do histograma de trend='cte' parece conter uma tendência, ### mas os modelos de 'retirada' usados deixam os dados sem uma estrutura de ### variação. Dificilmente vai conseguir ajustar um bom modelo nesses casos, pois ### o comportamento é muito errático. É possível que O erro que você estava obtendo ### seja decorrente da fragilidade do modelo ajustado. ### Realmente não ficou bom! Mas e aí? O que dá pra fazer? Tenta ssim mesmo? ### Ainda tem que testar a anisotropia, que não tem no seu script original. ?variog4 va <- variog4(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical") par(par.ori) plot(va) ### Pronto! Apareceu anisotropia... ### Tem duas sugestões de leitura abaixo pra dar uma ajuda... browseURL(' http://www.rc.unesp.br/igce/aplicada/DIDATICOS/LANDIM/tkrigagem.pdf') browseURL(' https://www.soils.org/publications/sssaj/abstracts/61/1/SS0610010298') ### Você vai ter que estudar adequadamente a questão da anisotropia em seus dados ### Isso não te isenta do efeito combinado anisotropia + tendência ### Só pra testar as possibilidades, vamos fazer o variograma na direção 135 graus ### Vai ter que entrar a direções, sendo mais fácil usar PI/n ### Tabela das direções (PI/n) para variogramas denom <- c(1:4,6,8,10,12,24,36,45,60); k <- (pi/180) ### Denominadores data.frame(rad=paste0('pi/', denom), rad2=round(pi/denom,4), deg=pi/denom/k) ### 135 é 3*45 ou 3*pi/4 va.135 <- variog(dados, uvec=uvec, pairs.min=pairs, dir=3*pi/4, max.dist=hmax, estimator="classical") plot(va.135) ### Aqui parece ser mais fácil ajustar. Talvez tenha que ajustar o 'lag' ### Avaliar a questão da distância entre amostras... ### Lembra do resíduo que separei anteriormente? ### Se fosse usar... va.r = variog(dados.res, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical") va.r.env <- variog.mc.env(dados.res, obj.var = va.r) plot(va.r, xlab="Dist. (m)", ylim=c(0,summary(va.r$v)[6]*1.1), env=va.r.env, main="Residual") ### Vou parar por aqui, porque antes de prosseguir vai precisar ver ### adequadamente essa questão da anisotropia. ### No caso do efeito ser só da anisotropia, vai simplificar a krigagem e você ### não deve mais ter o problema que vinha tendo. ### <END> Éder Comunello <c <comunello.eder(a)gmail.com>omunello.eder(a)gmail.com> Dourados, MS - [22 16.5'S, 54 49'W]

2 1

Recusão no R, problemas ao implementar o algoritimo insertionsort no R
by Augusto Ribas 29 Out '13

29 Out '13

Ola pessoa, eu estava tentando implementar o algorítimo insertion sort no R de forma recursiva, mas não estou conseguindo e também não consigo identificar onde estou errado. #Veja so. Eu tenho um vetor qualquer: vetor<-sample(100) vetor #Iterativamente ele funciona beleza, eu fiz assim: insertionsort<-function(vetor){ n<-length(vetor) for(i in 2:n) { aux=vetor[i] j=i-1 while(j>=1 && vetor[j]>aux) { vetor[j+1]<-vetor[j] j=j-1 } vetor[j+1]=aux } return(vetor) } insertionsort(vetor) #Mas quando tento usar recursão, não consigo fazer funcionar, ele não organiza insertionsortR<-function(vetor,n) { if(n>1) { insertionsortR(vetor,n-1) aux=vetor[n] i=n while(vetor[i-1]>aux && i>1) { vetor[i]<-vetor[i-1] i<-i-1 } vetor[i]<-aux } return(vetor) } insertionsortR(vetor,length(vetor)) #Usando o rcpp ele funcionou dessa forma aqui: cppFunction(" NumericVector insertionsortRC(NumericVector vetor, int n) { double aux; int i; if(n>1) { insertionsortRC(vetor,n-1); aux=vetor[n-1]; i=n-1; while(vetor[i-1]>aux && i>=0 ) { vetor[i]=vetor[i-1]; i--; } vetor[i]=aux; } return vetor; } ") insertionsortRC(vetor,length(vetor)) Alguém sabe dizer onde estou errado la no insertionsort no R. Eu imagino se estou imaginado errado o funcionamento dos vetores no R, se ele ta refazendo um novo vetor a cada chamada e no final nada da certo do jeito que estou fazendo. Recursão funciona beleza no R, por exemplo usando fatorial: fatorial<-function(n) { if(n==1) { return(1) } else { return(n*fatorial(n-1)) } } fatorial(5) Mas no caso do vetor, não entendi o funcionamento ali dos vetores, se alguém puder dar uma luz :) -- Grato Augusto C. A. Ribas Site Pessoal: http://recologia.com.br/ <http://augustoribas.heliohost.org> Github: https://github.com/Squiercg Lattes: http://lattes.cnpq.br/7355685961127056

2 2

Selecionar o primeiro numero negativo depois do último positivo.
by Alisson Lucrecio 29 Out '13

29 Out '13

Caro Colegas da lista r-br, Boa tarde. Eu estou precisando encontrar o primeiro numero negativo depois do último numero positivo de uma sequencia c(10, 12, 11, 5, 4, -5, -10, -11, -15, 4, 5, 6, 8, 10, 8, 5, 4, -6, -2, 4, 8). Assim o numero selecionado seriam -5 e -6. Algum de vocês poderia me ajudar com uma função que tivesse tal finalidade? Obrigado. Alisson Lucrécio da Costa

3 3

Não consigo fazer uma ANOVA com arquivo com 1e+06 registros.
by Luiz Roberto Martins Pinto 29 Out '13

29 Out '13

Caros companheiros da R-BR. Não consigo fazer uma ANOVA com arquivo com 1e+06 registros. Então... preciso de ajuda!!! Dados: http://www1.datafilehost.com/d/c0d31775 Meu pc R for windows 2.15.1(x64) 8 Gb de Memo load('RCBD_Data.Rdata') # Arquivo com 1e+06 registros m=aov(y~Treat+block,data=Data) summary(m) Mensagem de erro: Erro: não é possível alocar vetor de tamanho 14.9 Gb Além disso: Mensagens de aviso perdidas: 1: In model.matrix.default(mt, mf, contrasts) : Reached total allocation of 8086Mb: see help(memory.size) 2: In model.matrix.default(mt, mf, contrasts) : Reached total allocation of 8086Mb: see help(memory.size) 3: In model.matrix.default(mt, mf, contrasts) : Reached total allocation of 8086Mb: see help(memory.size) 4: In model.matrix.default(mt, mf, contrasts) : Reached total allocation of 8086Mb: see help(memory.size) Luiz Roberto Martins Pinto Prof. Pleno/DCET/UESC Laboratório de Estatística Computacional Universidade Estadual de Santa Cruz Ilhéus-Bahia luizroberto.uesc(a)gmail.com skype: lrmpinto http://lattes.cnpq.br/2732314327604831

6 24

Problemas com a atualização e instalação de pacotes no R no UBUNTU 13.04
by Luciano F. Sgarbi 29 Out '13

29 Out '13

Prezados, Sou novo no uso do OS linux e estou com problemas na atualização do R (versão 3.0.2) e na instalação de qualquer pacote no R. Meu sistema operacional é o UBUNTU 13.04, 32-bit. Se alguém puder dar uma força serei grato. Abraços, Luciano -- *Luciano F. Sgarbi * * Currículo Lattes <http://lattes.cnpq.br/9788198135883354> **Dep. Ecologia, Instituto de Ciências Biológicas I** * *Universidade Federal de Goiás, campus II, Goiânia-GO, Brasil*

2 5

Polinômios ortogonais
by Fernando Souza 28 Out '13

28 Out '13

Caros amigos, estou realizando uma regressão linear com contrastes polinomiais par obter o grau do polinômio que melhor se ajusta aos dados. Realizei o modelo no R e fiz também manualmente utilizando a tabela de divisores, multiplicadores para comparações ortogonais. Acontece que embora a análise de variância indique ser sigignificativo o modelo de 2ª ordem quando obtenho o modelo pelo método summary a regressão calculada manualmente difere daquela apresentada pelo R alguem pode me explicar o que está acontecendo. Segue o CMR: dados<-structure(list(Bloco = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L), .Label = c("1", "2", "3", "4", "5", "6"), class = "factor"), Espacamento = c(1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5), Producao = c(3.02, 2.795, 2.966, 2.552, 2.822, 3.334, 3.1, 2.651, 2.687, 2.543, 3.065, 2.741, 2.624, 2.84, 2.597, 3.109, 2.85, 2.759, 2.903, 2.57, 3.181, 2.759, 2.759, 2.525, 2.66, 3.244, 2.723, 2.507, 2.417, 3.002)), .Names = c("Bloco", "Espacamento", "Producao"), row.names = c(NA, -30L), class = "data.frame") dados$Bloco<-factor(dados$Bloco) dados$Espacamento<-factor(dados$Espacamento) modelo3<-lm(Producao~Bloco+ordered(Espacamento), data=dados) O modelo obtido manualmente foi: Y =4,721 - 0,0458X + 0,00025X² atenciosamente -- ===================================================== Fernando Souza Zootecnista, DSc. Produção Animal cel: (+55) 82 8113-8781 e-mail:nandodesouza@gmail.com ==================================================

3 3

[OFF] Mensagem duplicada!
by Andre Oliveira 27 Out '13

27 Out '13

Pessoal a lista se encontra com algum problema? Minhas duas ultimas postaguens foi enviada duas vezes na mesma mensagem sendo que enviei somente uma vez! André Oliveira Souza

1 0

Instalar pacotes!
by Andre Oliveira 27 Out '13

27 Out '13

Pessoas, preciso da ajuda de vocês estou tentado instalar estes pacotes e obtenho a aviso. package ‘divisors’ is not available (for R version 3.0.2) install.packages("descr"), dep=T) install.packages("itm", dep=T) install.packages("divisors",dep=T) Se eles não estão disponíveis para o R, minha versão consta a informação. andre@andre-Aspire-5920:~$ R -version ATEN??O: op??o desconhecida '-version' R version 3.0.2 (2013-09-25) -- "Frisbee Sailing" Copyright (C) 2013 The R Foundation for Statistical Computing Platform: i686-pc-linux-gnu (32-bit) R é um software livre e vem sem GARANTIA ALGUMA. Você pode redistribuí-lo sob certas circunstâncias. Digite 'license()' ou 'licence()' para detalhes de distribuição. R é um projeto colaborativo com muitos contribuidores. Digite 'contributors()' para obter mais informações e 'citation()' para saber como citar o R ou pacotes do R em publicações. Digite 'demo()' para demonstrações, 'help()' para o sistema on-line de ajuda, ou 'help.start()' para abrir o sistema de ajuda em HTML no seu navegador. Digite 'q()' para sair do R. [Área de trabalho anterior carregada] O que devo proceder para resolver? obrigado André Oliveira Souza

1 0

Ler banco de dados .xls
by Andre Oliveira 27 Out '13

27 Out '13

Boa tarde, eu sempre usei no linux os pacotes. 1) gdata 2) RODBC Existe algum outro pacote mais moderno e eficiente no momento? André Oliveira Souza

3 2