Abril 2014 - R-br - DINF UFPR

Problema com abordagem inflacionada de zeros
by ASANTOS 29 Abr '14

29 Abr '14

Boa dia Pessoal, Tenho alguns dados que estou com dificuldade de encontrar uma abordagem de análise, tenho uma variável n.insetos (são dados de contagem de insetos, mas tenho várias parcelas com 0 insetos) e a variável resposta n.ataques (são dados de contagem do número de plantas atacadas, também com zeros), o meu problema é que: quando não tenho insetos, não tenho ataque, mas também tenho as vezes insetos, mas não tenho ataque. Tentei pela abordagem de Poisson inflacionado por zeros, sem sucesso (minha variável explicativa tem inúmeros zeros e é uma informação com implicação biológica que não pode ser descartada). Alguém saberia me dizer que abordagem eu poderia tentar para modelar essa relação? Segue CRM com dados parciais e a tentativa de ajuste inflacionada por zeros: #------------------------------------------------------------------------------- # Definições da sessão. rm(list=ls()) require(pscl) require(VGAM) #------------------------------------------------------------------------------- #Dados n.insetos<-c(4,3,1,2,2,2,1,4,1,1,1,0,2,2,3,3,2,6,2,3,3,4,3,0,2,5,6,3,3,3,5,2,1, 2,1,1,0,2,2,1,3,6,3,1,2,0,3,1,0,4,2,2,3,2,0,2,2,2,1,2,4,4,3,5,2,1,5,2,3,2,5,3,3, 3,4,3,2,5,3,4,1,2,2,2,1,2,3,1,4,2,4,4,1,2,2,9,3,1,2,2,2,1,3,3,3,1,1,2,5,2,1,1,0, 5,0,4,0,3,2,2,3,4,1,1,2,1,1,2,3,3,2,5,2,2,2,3,1,2,6,1,2,3,1,2) n.ataques<-c(0,1,3,0,1,2,1,1,2,2,0,0,1,0,1,1,0,0,2,0,1,0,2,0,0,1,3,0,1,0,0,1,1, 1,0,1,0,0,1,1,3,1,1,1,2,0,1,1,0,0,0,1,0,0,0,1,1,0,0,1,2,0,1,2,0,1,1,2,0,1,0,0,1, 2,1,1,3,1,3,2,0,1,0,0,1,0,2,2,2,2,2,0,1,0,2,3,2,2,3,0,1,2,1,1,0,1,2,1,2,0,1,1,0, 1,0,0,0,1,0,0,1,0,1,0,0,2,2,3,1,0,2,1,4,1,1,0,0,3,0,2,0,2,1,0) plot(n.insetos, n.ataques) #------------------------------------------------------------------------------- #------------------------------------------------------------------------------- # Modelo completo. compl.mod <- zeroinfl(n.ataques~n.insetos|n.insetos) coef(compl.mod) summary(compl.mod) #------------------------------------------------------------------------------- Obrigado, -- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos(a)cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ======================================================================

2 1

for e If
by Evelyn Alecrim 28 Abr '14

28 Abr '14

Boa tarde, Tenho a seguinte situação, tenho três variáveis: Sit, d1 e d2. Gostaria que um valor (k) fosse calculado levando em consideração o valor de Sit, ou seja, para cada Sit uma fórmula diferente seria usada, como abaixo. Sit=c(3, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1) d1=c(1.27, 5.80, 2.42, 3.11, 5.25, 4.74, 5.42, 5.14, 5.27, 5.61, 6.09, 6.38, 6.34, 8.60, 7.36, 9.96, 10.86, 11.35, 8.94, 9.19, 8.66, 8.98, 9.86, 9.94, 11.59, 15.61, 6.28, 7.83, 7.94, 7.67, 8.62, 9.43, 9.00) d2=c(9.12, 10.81, 10.21, 7.11, 7.23, 7.62, 7.69, 7.59, 5.73, 5.39, 5.35, 5.42, 5.32, 3.84, 2.89, 0.70, 1.65, 2.32, 4.01, 4.16, 4.32, 4.74, 5.12, 3.01, 1.91, 7.33, 3.81, 2.27, 2.22, 4.44, 4.76, 7.35, 7.35) for(i in 1:seq_along(Sit)){ if(Sit==1){ k<-(100 + d1 ^ 2 - d2 ^ 2) / 20 } else if(Sit==2){ k<-(10 + (100 + d1 ^ 2 - d2 ^ 2) / 20) } else { k<-(400 + d1 ^ 2 - d2 ^ 2) / 40 } print(k) } No entanto, o R me retorna a seguinte mensagem: [1] 7.960963 7.919597 7.540307 8.978000 9.382240 9.110080 9.256007 9.220288 9.873500 10.060500 10.211640 10.283200 10.297330 11.480360 11.145437 [16] 12.467790 12.880427 13.086003 11.596087 11.678762 11.408330 11.454320 11.775130 12.243588 13.267000 14.748580 10.623057 11.403900 11.452880 10.977882 [31] 11.291170 10.872560 10.674437 Mensagens de aviso perdidas: 1: In 1:seq_along(Sit) : expressão numérica tem 33 elementos: somente o primeiro usado 2: In if (Sit == 1) { : a condição tem comprimento > 1 e somente o primeiro elemento será usado 3: In if (Sit == 2) { : a condição tem comprimento > 1 e somente o primeiro elemento será usado Parece que ele está calculando usando apenas a terceira fórmula. Alguém pode dar uma luz? Obrigada, Emily

7 8

Estimador de Máxima verossimilhança para parâmetros da dist. q-exponencial
by Romero Luiz M. Sales Filho 27 Abr '14

27 Abr '14

Caros amigos, estou precisando encontrar os estimadores para os parâmetros da dist. q-ponencial, cuja pdf é dada da seguinte forma: [image: { (2-q) \lambda e_q^{-\lambda x}}] onde, [image: e_q(x) = [1+(1-q)x]^{1 \over 1-q}]. Estou tendo problemas pois utilizando o pacote Optim do R, não ocorre convergência em alguns métodos, na verdade só há convergência usando o método SANN, e mesmo assim as estimativas estão muito diferentes dos parâmetros que foram utilizados para gerar uma amostra aleatória da dist. q-exp. O código utilizado no R é o seguinte: ## Gerando Valores de uma q-exponencial n<-10000 u<-runif(n) z<-1.8 beta<- 1.2 dados1 = (1/(beta*(z-2)))*((u^((1-z)/(2-z))-1)/((1-z)/(2-z))) ## Definindo a função de verossimilhança vero <- function(par,x){ q = par[1] lambda = par[2] ((sum(log(1 - ((1 - q)*x*lambda))))/(1 - q)) + ((n)*log(2 - q)) + (n*log(lambda)) } ## Otimizando os parâmetros - o método SANN sempre converge. optim(par=c(1.8,4),fn=vero, method="SANN",x=dados1 ) E a saída é a seguinte: $par [1] 1.000014 3.807168 $value [1] -284366598 $counts function gradient 10000 NA $convergence [1] 0 $message NULL Houve 50 ou mais avisos (use warnings() para ver os primeiros 50) Alguém pode me ajudar nessa questão, uma vez que as estimativas dos parâmetros estão muito diferentes dos parâmetros originais? Obrigado! Romero.

5 4

Erro: estouro da pilha de contexto (O que fazer?)
by Alexandre Loures 26 Abr '14

26 Abr '14

Boa noite! Estou usando a rotina abaixo para construir uma dummy para as fronteiras dos estados brasileiros. Porém, quando tendo rodar o R reporta o seguinte erro: Erro: estouro da pilha de contexto na linha 18 O que devo fazer? mydata$contiguity<-ifelse((origin=="ac") & (destination=="am"),1, ifelse((origin=="am") & (destination=="ro"),1,ifelse((origin=="am") & (destination=="rr"),1,ifelse((origin=="am") & (destination=="pa"),1,ifelse((origin=="am") & (destination=="mt"),1, ifelse((origin=="ro") & (destination=="mt"),1,ifelse((origin=="ro") & (destination=="am"),1, ifelse((origin=="rr") & (destination=="am"),1,ifelse((origin=="rr") & (destination=="pa"),1, ifelse((origin=="pa") & (destination=="mt"),1,ifelse((origin=="pa") & (destination=="am"),1,ifelse((origin=="pa") & (destination=="am"),1,ifelse((origin=="pa") & (destination=="rr"),1, ifelse((origin=="pa") & (destination=="ap"),1,ifelse((origin=="pa") & (destination=="ma"),1,ifelse((origin=="pa") & (destination=="to"),1,ifelse((origin=="pa") & (destination=="mt"),1, ifelse((origin=="ap") & (destination=="pa"),1, ifelse((origin=="ma") & (destination=="pa"),1,ifelse((origin=="ma") & (destination=="to"),1,ifelse((origin=="ma") & (destination=="pi"),1, ifelse((origin=="pi") & (destination=="ma"),1,ifelse((origin=="pi") & (destination=="to"),1,ifelse((origin=="pi") & (destination=="pe"),1,ifelse((origin=="pi") & (destination=="ce"),1, ifelse((origin=="ce") & (destination=="pi"),1,ifelse((origin=="ce") & (destination=="pe"),1,ifelse((origin=="ce") & (destination=="pb"),1,ifelse((origin=="ce") & (destination=="rn"),1, ifelse((origin=="rn") & (destination=="ce"),1,ifelse((origin=="rn") & (destination=="pb"),1, ifelse((origin=="pb") & (destination=="rn"),1,ifelse((origin=="pb") & (destination=="ce"),1,ifelse((origin=="pb") & (destination=="pe"),1, ifelse((origin=="pe") & (destination=="pb"),1,ifelse((origin=="pe") & (destination=="ce"),1,ifelse((origin=="pe") & (destination=="pi"),1,ifelse((origin=="pe") & (destination=="ba"),1, ifelse((origin=="pe") & (destination=="al"),1, ifelse((origin=="al") & (destination=="pe"),1,ifelse((origin=="al") & (destination=="ba"),1,ifelse((origin=="al") & (destination=="se"),1, ifelse((origin=="se") & (destination=="al"),1,ifelse((origin=="se") & (destination=="ba"),1, ifelse((origin=="ba") & (destination=="se"),1,ifelse((origin=="ba") & (destination=="al"),1,ifelse((origin=="ba") & (destination=="pe"),1,ifelse((origin=="ba") & (destination=="pi"),1, ifelse((origin=="ba") & (destination=="to"),1,ifelse((origin=="ba") & (destination=="go"),1,ifelse((origin=="ba") & (destination=="mg"),1,ifelse((origin=="ba") & (destination=="es"),1, ifelse((origin=="to") & (destination=="pa"),1,ifelse((origin=="to") & (destination=="mt"),1,ifelse((origin=="to") & (destination=="go"),1,ifelse((origin=="to") & (destination=="ba"),1, ifelse((origin=="to") & (destination=="pi"),1,ifelse((origin=="to") & (destination=="ma"),1, ifelse((origin=="mt") & (destination=="ro"),1,ifelse((origin=="mt") & (destination=="am"),1,ifelse((origin=="mt") & (destination=="pa"),1,ifelse((origin=="mt") & (destination=="to"),1, ifelse((origin=="mt") & (destination=="go"),1,ifelse((origin=="mt") & (destination=="ms"),1, ifelse((origin=="ms") & (destination=="mt"),1,ifelse((origin=="ms") & (destination=="go"),1,ifelse((origin=="ms") & (destination=="mg"),1,ifelse((origin=="ms") & (destination=="sp"),1, ifelse((origin=="ms") & (destination=="pr"),1, ifelse((origin=="go") & (destination=="mg"),1,ifelse((origin=="go") & (destination=="ba"),1,ifelse((origin=="go") & (destination=="to"),1,ifelse((origin=="go") & (destination=="mt"),1, ifelse((origin=="go") & (destination=="ms"),1,ifelse((origin=="go") & (destination=="df"),1, ifelse((origin=="mg") & (destination=="sp"),1,ifelse((origin=="mg") & (destination=="rj"),1,ifelse((origin=="mg") & (destination=="es"),1,ifelse((origin=="mg") & (destination=="ba"),1, ifelse((origin=="mg") & (destination=="go"),1,ifelse((origin=="mg") & (destination=="ms"),1, ifelse((origin=="sp") & (destination=="rj"),1,ifelse((origin=="sp") & (destination=="mg"),1,ifelse((origin=="sp") & (destination=="ms"),1,ifelse((origin=="sp") & (destination=="pr"),1, ifelse((origin=="rj") & (destination=="mg"),1,ifelse((origin=="rj") & (destination=="mg"),1,ifelse((origin=="rj") & (destination=="es"),1, ifelse((origin=="es") & (destination=="ba"),1,ifelse((origin=="es") & (destination=="mg"),1,ifelse((origin=="es") & (destination=="rj"),1, ifelse((origin=="pr") & (destination=="sp"),1,ifelse((origin=="pr") & (destintaion=="sc"),1,ifelse((origin=="pr") & (destination=="ms"),1, ifelse((origin=="sc") & (destination=="pr"),1,ifelse((origin=="sc") & (destination=="rs"),1, ifelse((origin=="rs") & (destination=="sc"),1, ifelse((origin=="df") & (destination=="go"),1,0)))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))) -- *Alexandre Rodrigues Loures*

2 2

Criar Dummy
by Alexandre Loures 25 Abr '14

25 Abr '14

Olá pessoal! Estou trabalhando numa base de dados e preciso criar algumas dummies, sendo que comando que geralmente utilizo é: mydata$dac<-ifelse(mydata$origin=="ac",1,0) Mas como preciso criar várias dummies, gostaria de saber se há como fazer isso de uma só vez? Isto é, utilizando um único comando. As dummies são para origin, destination e ano. Ou seja: Dummies para as regiões de origem (origin) Dummies para as regiões de destino (destination) Dummies para as regiões de origem com relação ao ano (Por exemplo: origin ac e ano 1998 = 1, ...) Dummies para as regiões de destino com relação ao ano (Por exemplo: destination ac e ano 1998 = 1, ...) Dummy para ano > str(mydata) 'data.frame': 4158 obs. of 9 variables: $ origin : Factor w/ 27 levels "ac","al","am",..: 1 1 1 1 1 1 1 1 1 1 ... $ destination: Factor w/ 78 levels "ac","africa",..: 3 3 5 5 7 7 13 13 17 17 ... $ ano : int 1998 1999 1999 1998 1999 1998 1999 1998 1998 1999 ... $ pib_i : num 1694329 1866542 1866542 1694329 1866542 ... $ pop_i : int 514050 527937 527937 514050 527937 514050 527937 514050 514050 527937 ... $ value : num 785 115 32290 0 0 ... $ pib_j : num 6676333 6935080 13933920 13198510 2128307 ... $ pop_j : num 2688117 2713203 2580860 2520684 439781 ... $ dist : num 5039 5039 1445 1445 NA ... Desde já muito obrigado! -- *Alexandre Rodrigues Loures*

3 5

Instalar pacotes no software R
by Jodavid Ferreira 25 Abr '14

25 Abr '14

Como posso instalar todos os pacotes de uma vez no software R no linux? -- Jodavid Ferreira Responsável pela comissão de Comunicação do MJ da RCCPB - http://rccpbjovem.com Membro da equipe Estadual do Ministério de Comunicação da RCCPB - http://rccpb.com.br Responsável pelo Ministério de Comunicação do Grupo de Oração Boa Esperança - http://grupoboaesperanca.com Bacharelando em Estatística - UFPB Estagiário no Laboratório de Estatística Aplicada ao Processamento de Imagens e Geoprocessamento - http://de.ufpb.br/~leapig Site Pessoal - http://jodavidferreira.com Twitter - @JodavidFerreira <http://twitter.com/JodavidFerreira>

5 6

[OFF-Topic] Distribuição Normal e teste de t
by Tropidurus Torquatus 25 Abr '14

25 Abr '14

Caros R-users, Coloquei como off-topic pois acredito que esta seja mais uma dúvida conceitual do que de R ppd. Dado que tenho dois conjuntos de dados: x <- c ( 13.772 , 10.119 , 11.224 , 21.999 , 17.181 , 27.477 , 21.693 , 26.706 , 14.576 , 19.689 , 12.494 , 19.511 , 16.785 , 16.551 , 20.780 , 0.000 , 25.867 , 28.223 , 26.059 , 21.668 , 18.209 , 20.186 , 28.320 , 26.946 , 16.218 , 16.814 , 16.859 , 0.000 , 0.000 , 3.986 , 1.625 , 17.286 , 28.149 , 19.404 , 21.669 , 22.058 , 14.947 , 14.512 , 21.177 , 29.997 , 20.958 , 15.123 , 28.932 , 23.807 , 30.589 , 10.208 , 20.076 , 18.897 , 14.666 , 13.955) y <- c ( 13.902 , 10.102 , 11.301 , 22.081 , 17.234 , 27.540 , 21.741 , 26.677 , 14.621 , 19.766 , 12.547 , 19.622 , 16.857 , 16.628 , 20.787 , 0.000 , 25.943 , 28.313 , 26.166 , 21.724 , 18.255 , 20.307 , 28.387 , 27.079 , 16.202 , 16.857 , 16.875 , 0.000 , 0.000 , 3.951 , 1.681 , 17.353 , 28.273 , 19.422 , 21.724 , 22.164 , 14.902 , 14.458 , 21.174 , 30.024 , 20.949 , 15.101 , 28.926 , 23.850 , 30.597 , 10.147 , 20.035 , 18.873 , 14.691 , 13.892) As duas amostras estão pareadas. As duas amostras falham para o teste de normalidade: > shapiro.test ( x ) Shapiro-Wilk normality test data: testetcompounds$tR...tM.UV.ABC W = 0.9373, p-value = 0.01057 > shapiro.test ( y ) Shapiro-Wilk normality test data: testetcompounds$tR...tM.UV.ABC.Group10 W = 0.9372, p-value = 0.01043 No livro Bayesian Computation with R Sugere que para distribuições não normais é necessário fazer o teste de Monte Carlo Com a seguinte função ======================================================= tstatistic = function ( y , x ) { m = length ( y ) n = length ( x ) sp = sqrt( ( ( m - 1 ) * sd ( y ) ^ 2 + ( n - 1 ) * sd ( x ) ^ 2 ) / ( m + n - 2 ) ) t = ( mean ( y ) - mean ( x ) ) / ( sp * sqrt ( 1 / m + 1 / n) ) return ( t ) } ======================================================== Retornando t = -0.02242564 Em seguida o teste: ======================================================== alpha = 0.01 ; m = 50 ; n = 50 # sets alpha, m and n ( alpha is the stated significance level, m and n are the samples sizes ) N = 10000 # sets the number of simulations n.reject = 0 # counter of num. of rejections for ( i in 1 : N ) { x = rnorm ( m , mean = x , sd = x ) # simulates xs from population 1 y = rnorm ( n , mean = y , sd = y ) # simulates ys from population 2 t = tstatistic ( y , x ) # computes the t statistic if ( abs ( t ) > qt ( 1 - alpha / 2 , n + m - 2) ) n.reject = n.reject + 1 # reject if |t| exceeds critical pt } true.sig.level = n.reject/N # est. is proportions of rejections ======================================================== Retornando true.sig.level = 0.0049 Problema Conceitual 1) Como interpretar esse resultado? Porém ainda há uma diferença entre o t calculado por essa fórmula e o t calculado com o comando t.stat (padrão do R) ================================================================ t = t.test ( y , x , mu = 0 , paired = TRUE , var.equal = TRUE , conf.level = 0.9999 ) t Paired t-test data: y and x t = 4.6437, df = 49, p-value = 2.594e-05 alternative hypothesis: true difference in means is not equal to 0 99.99 percent confidence interval: 0.003082235 0.067077765 sample estimates: mean of the differences 0.03508 ================================================================ E finalmente… minha média da diferencia é conhecida e não é zero, assim utilizei o seguinte comando: t = t.test ( y , x , mu = 0.87 , paired = TRUE , var.equal = TRUE , conf.level = 0.9999, alternative = "less” ) Paired t-test data: y and x t = -110.5231, df = 49, p-value < 2.2e-16 alternative hypothesis: true difference in means is less than 0.87 99.99 percent confidence interval: -Inf 0.06545182 sample estimates: mean of the differences 0.03508 2) Problema de Comando… como fazer a simulação de Monte Carlo na situação acima? * Obs.: aceito sugestões de leitura sobre o assunto e comando alternativos para comparar a médias Muito obrigado e me desculpe pelo longo e-mail. Thales

1 0

[Off] Pacotes/Artigos na Journal of Statistical Software
by FHRB Toledo 24 Abr '14

24 Abr '14

Recomendações do editor... interessante! http://www.r-bloggers.com/publishing-an-r-package-in-the-journal-of-statist… att, FH

1 0

gravando png
by Hélio Gallo Rocha 23 Abr '14

23 Abr '14

Caros, bom dia. Estou tentando aumentar a resolução de um plot utilizando o comando 'png'. Ele gera um arquivo, mas com 0kb e trava o visualizador de plot, ou seja, se tentar plotar outra imagem, ele não gera, tenho de fechar o trabalho Estou usando o comando png(filename = "k:/2014/imas/image.png", width = 1890, height = 2550, bg = "white", res = 300, restoreConsole = TRUE) Se alguém puder auxiliar, agradeço desde já -- Hélio Gallo Rocha IFSULDEMINAS - Câmpus Muzambinho

2 2

Help
by Bruce Kelly 23 Abr '14

23 Abr '14

Olá a todos, estou com um pequeno problema, estou usando o pacote 'rela' usando a função paf para aplica a técnica de analise fatorial, pois esta roda com NA, diferentemente da factanal. Bem quando tenho meu banco de dados .csv, tendo rodar mas ele da a seguinte mensagem de erro: Your dataset is not a numeric object, gostaria de tornar esse banco numérico, como posso fazer isso. Obrigado pela atenção de todos. Bruce Kelly

3 4