Março 2016 - R-br

SAR, CAR e GWR spatial models
by Wagner Tassinari 16 Mar '16

16 Mar '16

Olá pessoal, tudo bem ? Preciso trabalhar com os modelos espaciais SAR, CAR e GWR, para isso peguei um material do Luc Anselin para estudar a respeito do assunto. Observei que os modelos SAR e CAR podem ser ajustados por estes comandos: columbus.lag <- lagsarlm(CRIME ~ INC + HOVAL,data=columbus, col.listw) columbus.err <- errorsarlm(CRIME ~ INC + HOVAL,data=columbus, col.listw) columbus.gwr = gwr.sel(CRIME ~ INC + HOVAL, data=columbus, adapt=T) E posso verificar os valores preditos através dos comandos fitted.values Mas eu tenho uma amostra de 25 observações para poder validar o valor predito (Y_hat) destes modelos. Qual seria o comando para eu poder essas previsões ? Sera que existe um material que vocês indicariam para estudar ? Sei que no modelo GWR essa estimativa seria mais difícil, já que existirá um coef. de regressão para cada polígono. Muito obrigado, -- Wagner S. Tassinari Departamento de Matemática Universidade Federal Rural do Rio de Janeiro. BR-465, Km 7 - Seropedica, RJ - Brasil CEP: 23890-000 Skype: wagner.tassinari wtassinari(a)gmail.com tassinari(a)ufrrj.br ------------------------------------------------------- "Statistical thinking will one day be as necessary for efficient citizenship as the ability to read and write." (H.G.Wellis)

2 1

Poisson truncada
by Leonard Mendonça de Assis 15 Mar '16

15 Mar '16

Alguém sabe se já existe gerador de Poisson truncada no R? Preciso gerar uns números aleatórios segundo uma Poisson, mas não posso ter um número maior que X. Acredito que seja o caso de uma Poisson truncada. Leonard

3 2

Comportamento da Log-Verossimilhança - Dist. q-Exponencial
by Romero Luiz M. Sales Filho 15 Mar '16

15 Mar '16

Caros amigos, Estou trabalhando com a função de log-verossimilhança da dist. q-Exponencial. A partir dela desejo obter as estimativa dos parâmetros da distribuição. Utlizo o pacote optim do R juntamente com a função Nelder-Mead. Obtenho estimativas de uma amostra controlada, ou seja uma amostra que gero através do gerador de números aleatórios dessa distribuição e que portanto conheço os parâmetros que originaram a amostra. Para alguns valores dos Parâmetros a estimativa é bem próxima dos valores dos parâmetros, entretanto quando eu ploto o gráfico da função de log verossimilhança da distribuição que estou usando os pontos que maximizam os gráficos são diferentes daes estimativas encontradas. Gostaria de saber por que isso acontece??? Se eu otimizo a função através de máxima verossimilhança era esperado que o gráfico da função de log-verossimilhança apresentasse os pontos de máximo da função bem próximos das estimativas encontradas... Alguém pode me ajudar?? será que estou construindo a função de forma errada?? Já refiz tudo e creio que a função está escrita corretamente... O que pode está acontecendo... Abaixo escrevo o código em que que utilizo para gerar uma amostra aleatória, depois estimo os parâmetros dessa amostra gerada via Nelder-Mead (pacote optim) e por fim contruo o gráfico da verossimilhança... Note que os pontos que maximizam o gráfico diferem das estimativas de máxima verossimilhaça: ### Gerando a amostra baseando-se nos parâmetros iniciais ### Considerando um parâmetro q positivo. q=1.5 eta=5 n=1000 u<-runif(n) y = eta*((1-(u)^((1-q)/(2-q)))/(1-q)) ### Obtendo as estimativas dos parâmetros a partir da amostra gerada anteriormente vero <- function(par,x){ q = par[1] eta = par[2] saida<-((sum(log(1 - ((1 - q)*x*(1/eta)))))/(1 - q)) + ((n)*log(2 - q)) + (n*log(1/eta)) return(-saida) } saiday<-optim(par=c(1.8,4),fn=vero, method= "Nelder-Mead",x=y ) q_chap<-saiday[1]$par[1] eta_chap<-saiday[1]$par[2] q_chap eta_chap ### Construindo os gráficos de log-verossimilhança eta <-seq(1, 10, l = 25) q <- seq(1.1, 1.9, l = 25) f <- function(eta,q) { ((sum(log(1 - ((1 - q)*y*(1/eta)))))/(1 - q)) + (n*log(2 - q)) + (n*log(1/eta)) } z <- outer(eta, q, f) persp(eta, q, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue") Repito o mesmo código digitado acima, apenas com a alteração de que a amostra gerada agora é para um parâmetro q menor que 1, observe que a coisa fica ainda pior nesse caso, uma vez que, sempre o menor valor de q e o menor valor de eta são os pontos de máximo do gráfico: ### Gerando a amostra baseando-se nos parâmetros iniciais ### Considerando um parâmetro q menor que 1. q=0.5 eta=5 n=1000 u<-runif(n) y = eta*((1-(u)^((1-q)/(2-q)))/(1-q)) ### Obtendo as estimativas dos parâmetros a partir da amostra gerada anteriormente vero <- function(par,x){ q = par[1] eta = par[2] saida<-((sum(log(1 - ((1 - q)*x*(1/eta)))))/(1 - q)) + ((n)*log(2 - q)) + (n*log(1/eta)) return(-saida) } saiday<-optim(par=c(1.8,4),fn=vero, method= "Nelder-Mead",x=y ) q_chap<-saiday[1]$par[1] eta_chap<-saiday[1]$par[2] q_chap eta_chap ### Construindo os gráficos de log-verossimilhança eta <-seq(1, 10, l = 25) q <- seq(0, 0.75, l = 25) f <- function(eta,q) { ((sum(log(1 - ((1 - q)*y*(1/eta))),na.rm=TRUE))/(1 - q)) + (n*log(2 - q)) + (n*log(1/eta)) } z <- outer(eta, q, f) persp(eta, q, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue") Em tempo, coloco aqui algumas equações importantes: função PDF da dist. q-Exponencial: Verossimilhança da dist. q-Exponencial: log-verossimilhança da dist. q-Exponencial: gostaria muito de entender o que está acontecendo e por que nunca bate as estimativas dos parâmetros com o máximo do gráfico! Obrigado, Romero.

3 7

OFF TOPIC: 1o Congresso de Aplicações de Linguagem R em Administração COPPEAD/UFRJ.
by Marcos Silva 14 Mar '16

14 Mar '16

Talvez possa ser de interesse para alguns: http://cara.coppead.ufrj.br/ Abs. -- Marcos F. Silva http://sites.google.com/site/marcosfs2006

1 0

OFF-TOPIC: Aniversário da lista
by Éder Comunello 13 Mar '16

13 Mar '16

Pessoal, boa tarde! Em 15 de março a r-br completará 5 anos. Meus parabéns aos mantenedores da lista e todos seus usuários por essa importante fonte de aprendizado. Abaixo vai um código que encontrei na rede e que permite tocar o "Parabéns pra você" no R. ### <code r> # browseURL(" http://stackoverflow.com/questions/31782580/how-can-i-play-birthday-music-u… ") library("dplyr") library("audio") notes <- c(A = 0, B = 2, C = 3, D = 5, E = 7, F = 8, G = 10) pitch <- "D D E D G F# D D E D A G D D D5 B G F# E C5 C5 B G A G" duration <- c(rep(c(0.75, 0.25, 1, 1, 1, 2), 2), 0.75, 0.25, 1, 1, 1, 1, 1, 0.75, 0.25, 1, 1, 1, 2) bday <- data_frame(pitch = strsplit(pitch, " ")[[1]], duration = duration) bday <- bday %>% mutate(octave = substring(pitch, nchar(pitch)) %>% {suppressWarnings(as.numeric(.))} %>% ifelse(is.na(.), 4, .), note = notes[substr(pitch, 1, 1)], note = note + grepl("#", pitch) - grepl("b", pitch) + octave * 12 + 12 * (note < 3), freq = 2 ^ ((note - 60) / 12) * 440) tempo <- 120 sample_rate <- 44100 make_sine <- function(freq, duration) { wave <- sin(seq(0, duration / tempo * 60, 1 / sample_rate) * freq * 2 * pi) fade <- seq(0, 1, 50 / sample_rate) wave * c(fade, rep(1, length(wave) - 2 * length(fade)), rev(fade)) } bday_wave <- mapply(make_sine, bday$freq, bday$duration) %>% do.call("c", .) play(bday_wave) ### </code> ================================================ Éder Comunello Agronomist (UEM), MSc in Environ. Sciences (UEM) DSc in Agricultural Systems Engineering (USP/Esalq) Brazilian Agricultural Research Corporation (Embrapa) Dourados, MS, Brazil |<O>| ================================================ GEO, -22.2752, -54.8182, 408m UTC-04:00 / DST: UTC-03:00

1 0

GlobalEnv
by Wagner Tassinari 11 Mar '16

11 Mar '16

Olá pessoal, Sei que é besteira, mas o que significa esta mensagem abaixo, pois não consigo atachar o meu banco quando ela aparece. Como devo proceder ? > attach(data) The following object is masked _by_ .GlobalEnv: cadastro Obrigado pessoal... abs -- Wagner S. Tassinari Departamento de Matemática Universidade Federal Rural do Rio de Janeiro. BR-465, Km 7 - Seropedica, RJ - Brasil CEP: 23890-000 Skype: wagner.tassinari wtassinari(a)gmail.com tassinari(a)ufrrj.br ------------------------------------------------------- "Statistical thinking will one day be as necessary for efficient citizenship as the ability to read and write." (H.G.Wellis)

2 1

Criar comparações múltiplas entre fatores
by ASANTOS 11 Mar '16

11 Mar '16

Caros listeiros, Gostaria de achar uma maneira de automatizar a criação de comparações múltiplas entre fatores, no meu caso: #Fator treat <- gl(3, 20, labels = paste("t", 1:3, sep="")) #Variables sp1_1 = rnorm(10, 5, 0.25) sp1_2 = rnorm(50, 2.5, 0.25) sp1 = c(sp1_1,sp1_2) sp2 = rnorm(60, 2.5, 0.25) sp3_1 = rnorm(10, 12, 0.25) sp3_2 = rnorm(50, 2.5, 0.25) sp3 = c(sp3_1,sp3_2) sp4 = rnorm(60, 2.5, 0.25) sp<-cbind(sp1,sp2,sp3,sp4) #Criação da matriz de contrastes Treat_Imp<-model.matrix(~ treat-1) Treat_Imp1_2<-Treat_Imp[, 1] - Treat_Imp[, 2] Treat_Imp2_3<-Treat_Imp[, 2] - Treat_Imp[, 3] Treat_Imp1_3<-Treat_Imp[, 1] - Treat_Imp[, 3] Eu acabei criando manualmente as comparações, através dos objetos Treat_Imp1_2, Treat_Imp2_3 e Treat_Imp1_3. Mas alguém poderia me dar uma ideia de como sistematizar isso, e ter a criação desses objetos de uma maneira direta e que não envolva pacotes como o multcomp? 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 OrcID: orcid.org/0000-0001-8232-6722 Researchgate: https://www.researchgate.net/profile/Alexandre_Santos10 LinkedIn: https://br.linkedin.com/in/alexandre-dos-santos-87961635 ======================================================================

2 3

I Congresso de Aplicações de R em Administração
by Celina Rebello 11 Mar '16

11 Mar '16

O I Congresso de Aplicações de R em Administração - CARA - visa disseminar o universo de estudos na área de Administração que utilizam a linguagem R como forma de alavancar a eficiência e qualidade dos trabalhos elaborados. Haverá espaço para exposição de trabalhos das diversas áreas de Administração, que tenham sido desenvolvidos utilizando a linguagem R. O foco não deverá ser a discussão do código desenvolvido, mas sim a análise proposta, e do poder estatístico oferecido pela linguagem. Interessados devem submeter até três artigos pelo sistema disponível neste site, observando o calendário do evento. O CARA é destinado tanto a profissionais que nunca escreveram um código em R, bem como a profissionais experientes que desejarem conhecer trabalhos que têm sido desenvolvidos utilizando-se desta poderosa linguagem. Serão também oferecidas duas atividades: Workshop Análise de Business Intelligence, para profissionais que possuem um conhecimento básico de R; Minicurso Apresentação do R, voltado àqueles que desejarem ter o primeiro contato com a linguagem. Submeta seus trabalhos em: http://cara.coppead.ufrj.br/submissao-de-trabalhos Saiba Mais em :http://cara.coppead.ufrj.br/home

1 0

Bons valores iniciais e não convergência - regressão não linear
by Adriele Pereira 08 Mar '16

08 Mar '16

Olá, bom dia a todos. Estou tentando ajustar o modelo não linear Richards ao seguinte conjunto de dados: y<-c(32.2,45.1,63.3,87.7,107.6,120.4,131.8,148.9,159.3,166.8,173.4,187.4,195.5,198.8,211.0,221.4,226.6,231.2,237.7,243.1,246.6,250.3) x<-c(90,180,270,360,450,540,630,720,810,900,990,1080,1170,1260,1350,1440,1530,1620,1710,1800,1890,1980) Ao utilizar o shiny, obtive os seguintes valores iniciais: alpha=270 , beta= 402, k=0.0016 , d=0.8. O problema é que o gráfico fica muito bom com esses valores, mas ao colocar esse valores iniciais na nls aparece a mensagem de erro: ajuste=nls(y~(alpha*(1 + (d - 1)*exp(-k*(x-beta)))^(1/(1-d))), start=c(alpha= 270, beta= 402 , k= 0.0016, d=0.8)) Error in numericDeriv(form[[3L]], names(ind), env) : Obtido valor faltante ou infinito quando avaliando o modelo Já tentei mudar a parametrização do modelo e aumentar o número de iterações e mesmo assim não ajustou. Alguém já se deparou com um problema parecido? Desde já agradeço.

2 1

Soma com intervalo definido dentro de um loop
by Yury Duarte 07 Mar '16

07 Mar '16

Bom dia colegas de R! Preciso somar uma serie de dados com tamanho definido igual a "a" repetições. Tentei usar o comando for() para realizar essa tarefa e somar os valores de x1 até xa; x2 até x(a+1) e assim sucessivamente até xn. Usando a estrutura básica da função soma fora do comando for não tive problemas: sum(Nasa$PP.Liquida[1:140]) Quando tentei aplicar esse formato no comando for não consegui ter sucesso: PP = for(i in Nasa$PP.Liquida){ sum(Nasa$PP.Liquida[i:(i+139)]) } Acredito que possa ser um problema de entendimento conceitual do funcionamento desse comando. Alguém poderia me ajudar nessa questão? Desde já, agradeço a todos pelo apoio! Att Yury Duarte Engenheiro Agrônomo - ESALQ/USP

3 3