- R-br - DINF UFPR

Diferença na interpretação da diferença da log-verossimilhança e do teste de Vuong em modelos inflacionados de zeros
by ASANTOS 25 Jan '15

25 Jan '15

Caros Membros, Estou com uma dúvida referente a diferença na interpretação dos resultados da diferença da log-verossimilhança e do teste de Vuong em modelos inflacionados de zeros, ajustados com o pacote pscl. Fiz a junção de dois tratamentos e ajustei um novo modelo e comparei com o modelo completo pela diferença entre a log-verossimilhança dos modelos, sendo: > pchisq(D, df=df, lower.tail=FALSE)#Os tratamento 1 e 2 são diferentes [1] 4.986542e-07 e em seguida realizei a comparação pelo teste de Vuong (que utiliza a razão de verossimilhança entre modelos), sendo: Vuong z-statistic H_A p-value Raw 1.0541254 model1 > model2 0.14591 AIC-corrected 1.0517040 model1 > model2 0.14647 BIC-corrected 0.7567064 model1 > model2 0.22461 porém os resultados foram diversos e fiquei em dúvida. Como os tratamentos podem ser diferentes e ao mesmo tempo o modelo que os considera semelhantes é significativamente igual ao modelo completo? Segue CRM com meu problema com dados artificiais, obrigado, #------------------------------------------------------------------ # Definições da sessão rm(list=ls()) require(pscl) #------------------------------------------------------------------ # Dados artificiais da <- expand.grid(trat=gl(4,1), tempo=1:15) X <- model.matrix(~trat+tempo, da); ncol(X) betas <- c(0.1,0.3,0.5,0.7,0.8) eta <- X%*%betas y1 <- rpois(da$trat, lambda=exp(eta)) y2 <- rbinom(y1, size=1, prob=0.20) da$y <- y1*y2 str(da) #------------------------------------------------------------------ # Ajuste do modelo completo m1 <- zeroinfl(y~trat+tempo|trat, data=da) summary(m1) # sort(tapply(da$y,da$trat,mean))# Ordena as medias por tratamento # ########### Comparando pela junção dos tratamentos 1 e 2 ----------- trat_1_2<-da$trat levels(trat_1_2) levels(trat_1_2)[1]<-"trat_1_2" levels(trat_1_2)[2]<-"trat_1_2" levels(trat_1_2) # m2 <- zeroinfl(y~trat_1_2+tempo|trat_1_2, data=da) summary(m2) # # Comparando o modelo completo e o ajunte com junção dos tratamentos- df <- length(coef(m1))-length(coef(m2)) D <- 2*abs(diff(c(logLik(m1), logLik(m2)))) D pchisq(D, df=df, lower.tail=FALSE)#Os tratamento 1 e 2 são diferentes # # # Comparando os modelos pelo teste de Vuong------------------------- vuong(m1,m2)#Porém os modelos são iguais # #<END>---------------------------------------------------------------- -- ====================================================================== 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 ====================================================================== --- Este email foi escaneado pelo Avast antivírus. http://www.avast.com

1 0

GLM de Poisson para dados com ponto máximo
by ASANTOS 25 Jan '15

25 Jan '15

Bom dia Pessoal, Estou com problemas no ajuste de um modelo GLM de Poisson, trata-se de uma relação entre número de insetos (n.insetos) com a temperatura (temp). O ajuste usando o modelo GLM de Poisson (mod.1) não parece ser o mais adequado e demonstra sobredispersão, além que a interpretação da curva do ponto de vista biológico não é realistica, pois eu teria que ter uma curva que aumenta com a temperatura, chega a um ponto ótimo e diminui, como uma equação de 2º grau, porém o modelo de Poisson que utilizei não faz isto, por isso pergunto, existe alguma sugestão de como eu poderia ajustar dados dessa natureza com mais semelhança biológica com o fenômeno observado? Segue CRM como exemplo: ##Dados artificias # Número de indivíduos n.insetos<-c(rpois(10,lambda=5),rpois(10,lambda=8),rpois(10,lambda=20),rpois(10,lambda=30),rpois(10,lambda=7),rpois(10,lambda=5)) #Temperatura temp<- c(rnorm(10,10,0.1),rnorm(10,15,0.1),rnorm(10,22,0.1),rnorm(10,24,0.1),rnorm(10,28,0.1),rnorm(10,32,0.1)) # #Visualiza o comportamentos das variáveis plot(temp,n.insetos) #Ajuste por GLM de Poisson mod.1<-glm(n.insetos~temp, family="poisson") summary(mod.1) anova(mod.1, test="Chi")## Anova #Representação do modelo a<-mod.1$coeff[1]###Coeficiente a b<-mod.1$coeff[2] ###Coeficiente b x<-seq(5,38,1) curve(exp(a+b*x),lty=1, add=T) # 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 ====================================================================== --- Este email foi escaneado pelo Avast antivírus. http://www.avast.com

2 2

Microsoft adquire Revolution Analytics
by David Feitosa 23 Jan '15

23 Jan '15

Microsoft to acquire Revolution Analytics to help customers find big data value with advanced statistical analysis http://blogs.microsoft.com/blog/2015/01/23/microsoft-acquire-revolution-ana… Revolution Analytics joins Microsoft http://blog.revolutionanalytics.com/2015/01/revolution-acquired.html Atenciosamente, David Feitosa (\_(\ (=°;°) (("")("")

1 0

Ordenação de importância de variáveis dependentes usando permanova
by Alexandre Santos 23 Jan '15

23 Jan '15

Caros listeiros, Estou tentando ordenar as variáveis dependentes mais importantes na projeção, dentre de 5 (y1 a y5), responsáveis pelas diferenças apresentadas por três tratamentos (colony1, colony2 e colony3) usando Permanova com a função adonis do pacote vegan. Porém, não estou conseguindo extrair essa informação do objeto man.col2, criado pela função adonis para ordenar as variáveis dependentes mais importantes que no meu caso separaram os tratamentos colony 1 e 2 da colony 3, alguém com familiaridade com a função adonis poderia me ajudar, segue CRM: #Pacotesrequire(vegan) #Dados artificiais#Tratamentoscolony<-as.factor(sort(rep(c("colony1", "colony2","colony3"), 100)))#Variáveis dependentesy1 <- c(rnorm(100,1,0.1),rnorm(100,1,0.1),rnorm(100,3.5,0.1))y2 <- c(rnorm(100,50,0.3),rnorm(100,50,0.3),rnorm(100,5,0.6))y3 <- c(rnorm(100,10,2.3),rnorm(100,10,2.3),rnorm(100,11,2.5))y4 <- c(rnorm(100,5,0.5),rnorm(100,5,0.5),rnorm(100,22,0.5))y5 <- c(rnorm(100,15,0.1),rnorm(100,15,0.1),rnorm(100,14,0.1)) #Cria o data frameavaliacao <- as.factor(colony)espectro <- cbind(y1,y2,y3,y4,y5)dados <- data.frame(avaliacao = I(as.matrix(avaliacao)), bands = I(as.matrix(espectro))) ## Cria uma matriz de contraste contrasts(avaliacao) <- cbind(c(1, -1, 0), c(1, 0, -1))avaliacao.ortho <- model.matrix(~ avaliacao) ## Cria um fator para cada contrastecol1vs2 <-avaliacao.ortho[, 2]col1vs3<-avaliacao.ortho[, 3] #Permanova com adonis para o modelo completoman.col<-adonis(bands ~ avaliacao, data=dados,permutations=99)man.col #Permananova com adonis para o modelo com o modelo decomposto#Inserindo os contrastes de interesseman.col2 <- adonis(dados$bands ~ col1vs2 + col1vs3, permutations=99)man.col2## Obrigado-- ======================================================================Alexandre dos SantosProteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato GrossoCampus CáceresCaixa Postal 244Avenida dos Ramires, s/nBairro: Distrito Industrial Cáceres - MT CEP: 78.200-000Fone: (+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 =====================================================================

1 0

Ordenação de importância de variáveis dependentes usando permanova
by Alexandre Santos 23 Jan '15

23 Jan '15

Caros listeiros, Estou tentando ordenar as variáveis dependentes mais importantes na projeção, dentre de 5 (y1 a y5), responsáveis pelas diferenças apresentadas por três tratamentos (colony1, colony2 e colony3) usando Permanova com a função adonis do pacote vegan. Porém, não estou conseguindo extrair essa informação do objeto man.col2, criado pela função adonis para ordenar as variáveis dependentes mais importantes que no meu caso separaram os tratamentos colony 1 e 2 da colony 3, alguém com familiaridade com a função adonis poderia me ajudar, segue CRM: #Pacotesrequire(vegan) #Dados artificiais#Tratamentoscolony<-as.factor(sort(rep(c("colony1", "colony2","colony3"), 100)))#Variáveis dependentesy1 <- c(rnorm(100,1,0.1),rnorm(100,1,0.1),rnorm(100,3.5,0.1))y2 <- c(rnorm(100,50,0.3),rnorm(100,50,0.3),rnorm(100,5,0.6))y3 <- c(rnorm(100,10,2.3),rnorm(100,10,2.3),rnorm(100,11,2.5))y4 <- c(rnorm(100,5,0.5),rnorm(100,5,0.5),rnorm(100,22,0.5))y5 <- c(rnorm(100,15,0.1),rnorm(100,15,0.1),rnorm(100,14,0.1)) #Cria o data frameavaliacao <- as.factor(colony)espectro <- cbind(y1,y2,y3,y4,y5)dados <- data.frame(avaliacao = I(as.matrix(avaliacao)), bands = I(as.matrix(espectro))) ## Cria uma matriz de contraste contrasts(avaliacao) <- cbind(c(1, -1, 0), c(1, 0, -1))avaliacao.ortho <- model.matrix(~ avaliacao) ## Cria um fator para cada contrastecol1vs2 <-avaliacao.ortho[, 2]col1vs3<-avaliacao.ortho[, 3] #Permanova com adonis para o modelo completoman.col<-adonis(bands ~ avaliacao, data=dados,permutations=99)man.col #Permananova com adonis para o modelo com o modelo decomposto#Inserindo os contrastes de interesseman.col2 <- adonis(dados$bands ~ col1vs2 + col1vs3, permutations=99)man.col2## Obrigado-- ======================================================================Alexandre dos SantosProteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato GrossoCampus CáceresCaixa Postal 244Avenida dos Ramires, s/nBairro: Distrito Industrial Cáceres - MT CEP: 78.200-000Fone: (+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 =====================================================================

1 0

Comparação de slope e intercepto entre regressões lineares
by marcelo claro de souza 22 Jan '15

22 Jan '15

Bom dia pessoal Eu plotei 10 regressões lineares, tipo dose resposta, e gostaria de comparar se elas diferem entre si em relação ao slope ou intercepto. Alguma dica de script para fazer uma comparação múltipla entre as regressões? Muito obrigado pela colaboração. Abraço Marcelo -- Marcelo Claro de Souza Post doc at the University of São Paulo School of Pharmaceutical Sciences of Ribeirão Preto Laboratory of Pharmacognosy Brazil

2 2

Ajuda com a criação de um looping
by Simone D. Sartorio 21 Jan '15

21 Jan '15

Pessoal, preciso de uma ajuda... Meus dados são esses "a21", nele tenho a observação de 219 vacas (n_vaca), que estão numeradas de maneira aleatória. As vacas foram separadas em 3 grupos, onde em cada grupo aplicaram um tratamento diferente. Esse tratamento foi aplicado em forma de injeções (17 injeções) e entre uma injeção e outra, temos 14 dias de intervalo. A variável de interesse é "producao"> head(a21) n_vaca n_lactacao trat pri_multiparas n_injecao dia_ciclo dia_exp producao saude8 1731 3 1 2 0 1 -13 45.1 healthy9 1731 3 1 2 0 2 -12 47.8 healthy10 1731 3 1 2 0 3 -11 45.5 healthy11 1731 3 1 2 0 4 -10 45.2 healthy12 1731 3 1 2 0 5 -9 47.4 enrollment13 1731 3 1 2 0 6 -8 48.3 enrollment # Aqui estou criando um subconjunto, onde tenho todas as informações da vaca 1731:> (vaca1731<- subset(a21, a21$n_vaca==1731)); dim(vaca1731) n_vaca n_lactacao trat pri_multiparas n_injecao dia_ciclo dia_exp producao saude8 1731 3 1 2 0 1 -13 45.1 healthy9 1731 3 1 2 0 2 -12 47.8 healthy10 1731 3 1 2 0 3 -11 45.5 healthy11 1731 3 1 2 0 4 -10 45.2 healthy12 1731 3 1 2 0 5 -9 47.4 enrollment13 1731 3 1 2 0 6 -8 48.3 enrollment14 1731 3 1 2 0 7 -7 46.4 enrollment15 1731 3 1 2 0 8 -6 47.9 enrollment16 1731 3 1 2 0 9 -5 50.3 enrollment17 1731 3 1 2 0 10 -4 50.8 enrollment18 1731 3 1 2 0 11 -3 47.9 enrollment19 1731 3 1 2 0 12 -2 46.2 heat20 1731 3 1 2 0 13 -1 50.8 healthy21 1731 3 1 2 0 14 0 45.8 healthy22 1731 3 1 2 1 1 1 50.5 healthy23 1731 3 1 2 1 2 2 49.0 healthy24 1731 3 1 2 1 3 3 48.8 healthy25 1731 3 1 2 1 4 4 53.9 healthy26 1731 3 1 2 1 5 5 45.4 healthy27 1731 3 1 2 1 6 6 38.6 healthy28 1731 3 1 2 1 7 7 41.3 healthy29 1731 3 1 2 1 8 8 50.2 healthy30 1731 3 1 2 1 9 9 50.2 healthy31 1731 3 1 2 1 10 10 50.4 healthy32 1731 3 1 2 1 11 11 49.9 healthy33 1731 3 1 2 1 12 12 50.1 healthy34 1731 3 1 2 1 13 13 52.0 healthy35 1731 3 1 2 1 14 14 47.8 healthy36 1731 3 1 2 2 1 15 50.0 healthy37 1731 3 1 2 2 2 16 49.7 healthy38 1731 3 1 2 2 3 17 50.4 healthy39 1731 3 1 2 2 4 18 48.9 healthy40 1731 3 1 2 2 5 19 49.5 healthy41 1731 3 1 2 2 6 20 48.8 healthy42 1731 3 1 2 2 7 21 56.3 healthy43 1731 3 1 2 2 8 22 52.1 healthy44 1731 3 1 2 2 9 23 51.6 healthy45 1731 3 1 2 2 10 24 50.5 healthy46 1731 3 1 2 2 11 25 46.0 healthy47 1731 3 1 2 2 12 26 53.7 healthy48 1731 3 1 2 2 13 27 51.7 healthy49 1731 3 1 2 2 14 28 49.1 healthy50 1731 3 1 2 3 1 29 51.9 healthy51 1731 3 1 2 3 2 30 51.4 healthy...256 1731 3 1 2 17 11 235 21.2 lamness257 1731 3 1 2 17 12 236 19.4 lamness258 1731 3 1 2 17 13 237 23.5 lamness259 1731 3 1 2 17 14 238 19.9 lamness[1] 251 9 # Para simplificar, considere que temos apenas 10 vacas:#--------------------------------------------> vacas <- c(1731, 2090, 2652, 2945, 2946, 3246, 3413, 3593, 3878, 3882) # é a numeração das vacas # O que eu quero fazer com esses dados é o seguinte:v1<- subset(a21, a21$n_vaca==vacas[1]); dim(v1) v2<- subset(a21, a21$n_vaca==vacas[2]); dim(v2)v3<- subset(a21, a21$n_vaca==vacas[3]); dim(v3)v4<- subset(a21, a21$n_vaca==vacas[4]); dim(v4)v5<- subset(a21, a21$n_vaca==vacas[5]); dim(v5)v6<- subset(a21, a21$n_vaca==vacas[6]); dim(v6)v7<- subset(a21, a21$n_vaca==vacas[7]); dim(v7)v8<- subset(a21, a21$n_vaca==vacas[8]); dim(v8)v9<- subset(a21, a21$n_vaca==vacas[9]); dim(v9)v10<- subset(a21, a21$n_vaca==vacas[10]); dim(v10) b<- c(tapply(v1$producao, v1$dia_ciclo, mean),tapply(v2$producao, v2$dia_ciclo, mean),tapply(v3$producao, v3$dia_ciclo, mean),tapply(v4$producao, v4$dia_ciclo, mean),tapply(v5$producao, v5$dia_ciclo, mean),tapply(v6$producao, v6$dia_ciclo, mean),tapply(v7$producao, v7$dia_ciclo, mean),tapply(v8$producao, v8$dia_ciclo, mean),tapply(v9$producao, v9$dia_ciclo, mean),tapply(v10$producao, v10$dia_ciclo, mean))#-------------------------------------------- Ainda estou engatinhando em questão de programação no R, como poderia fazer um looping com essas linhas de comando, pois preciso fazer isso para as 219 vacas! desde já, muito obrigada!Simone ***---***---***---***---***---***---***---***---***---***---***---***---**** Simone Daniela Sartorio* Professora Adjunta I no CCA/UFSCar, Araras/SP.* Doutora e Mestre em Estatística e Experimentação Agronômica - ESALQ/USP.* Licenciada em Matemática - UNESP/Rio Claro.***---***---***---***---***---***---***---***---***---***---***---***---*** Tenha um bom dia! ;)

3 7

Pacote plyr - group_by
by Henry Mazer 21 Jan '15

21 Jan '15

Olá! No csv http://www.datafilehost.com/d/32b663ce tenho 4 variáveis.Quero saber quantas pessoas são homens, víuvos, com idade acima de 79 e se formaram no Paraná e também todos os outros cruzamentos. Fiz o seguinte:teste <- read.csv('cruzamento1.csv')library("dplyr")group_by(teste, Qual.o.seu.gênero., Qual.o.seu.estado.civil., grupos.de.idade, Em.qual.estado.você.se.formou.)%>%summarise(sum) E tive como resposta:Erro: not a vector O que fiz de errado? Agradeço desde já.

4 4

Gerar um executável de um arquivo.r
by Thaís Brenda Martins 21 Jan '15

21 Jan '15

Boa tarde pessoal, Meu trabalho de conclusão de curso foi desenvolver uma ferramenta gráfica para teste de hipótese sob médias via R. Mas estou com três scripts, e queria deixar com apenas um executável, ou seja, para que o usuário só clicasse no ícone e já executasse a interface gráfica. Atualmente eu tenho que passar os três arquivos.r pra pessoa e através do R ele consegue executar. Como posso fazer para compilar estes arquivos e gerar um executável? Att.,Thaís Brenda Martins

3 2

off-topic Criar função no emacs para fazer um Text-template
by Augusto Ribas 20 Jan '15

20 Jan '15

Ola pessoal. Estou tentando aprender como faz para criar funções no Emacs, mais especificamente macro. Algo do tipo, eu aperto ctrl c mais alguma coisa e ele colocar um trechinho de texto no buffer que ta selecionando. Por exemplo algo assim ##################### # ##################### Pra mim escrever um titulo, sei la. Alguém tem algum exemplo de algo assim? Eu estou apanhando um pouco de alguns manuais que estou tentando ler. Bem agradeço a atenção. -- 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

4 4