- R-br - DINF UFPR

Ajuda: HPD em conjunto de valores preditos
by Flavia Cristina M.Queiroz 22 Mai '13

22 Mai '13

Elias... obrigada por responder... Então não tem como eu calcular o HPD para meus dados? Acontece que eu estou achando a distribuição empírica das minhas predições e gostaria de obter um certo limite que me dê uma alta cobertura (um limite onde tenho máxima densidade)... Como eu poderia fazer isso?? Att. Flávia Cristina ---------------------------------------------------------------------- Message: 5 Date: Tue, 21 May 2013 18:21:57 -0300 From: Elias Teixeira Krainski <eliaskrainski(a)yahoo.com.br> To: r-br(a)listas.c3sl.ufpr.br Subject: Re: [R-br] Ajuda: HPD em conjunto de valores preditos Message-ID: <519BE575.5040900(a)yahoo.com.br> Content-Type: text/plain; charset="iso-8859-1"; Format="flowed" Faz sentido calcular HPD de uma distribuição, não de um grid de valores como vc fez. A função HPDinterval() aplica-se especificamente a amostras de uma distribuição. Você aplicou-a a uma sequencia de valores para os quais foi avaliada uma densidade. Ou seja, exatamente o mesmo que fazer: HPDinterval(as.mcmc(seq(-1, 1, length=512))) On 21/05/13 14:04, Flavia Cristina M.Queiroz wrote: > Olá pessoal, > Estou precisando de uma ajuda com o cálculo de intervalo de credibilidade: > > Tenho um conjunto de dados (valores preditos) com os quais tenho que > calcular o intervalo de credibilidade de predição. > Estes valores preditos foram calculados em uma implementação de redes > neurais. > > Daí eu calculei o HPD usando a função HPDinterval(), do pacote coda. > Fiz assim: > despad <- read.table('vintemaio_pad.txt',h=F) > dados <- as.matrix(despad) > a=density(dados,bw=0.45,from = -1, to = 1) > library(coda) > arq<-as.mcmc(a$x) > mcmcsamp(a$x) > HPD<-HPDinterval(arq) > HPD > Contudo, meu intervalo é muito amplo.. tipo, meus dados vão de -1 a 1 > e meu HPD deu [-0.9882583, 0.9138943]. > Porque meu intervalo deu tão amplo?? O que teria que fazer de diferente?? > ------------------------------------------------------------------------------- > Aguardo alguma dica ou solução!! > Att., > Flávia Cristina

2 1

Legenda externa do gráfico no base
by Alisson Lucrecio 22 Mai '13

22 Mai '13

Caro Colegas, Bom dia. Eu estou construindo 4 PCAs no R e gostaria de plota-las na mesma janela "par(mfrow=c(2,2))", entretanto não sei como fazer "par" ou mesmo redimensionar a janela para colocar a legenda no centro e no lado direito. Alguém sabe como fazer? Obrigado. Alisson Lucrécio da Costa

2 2

Classificar fosforo e textura do solo
by Hélio Gallo Rocha 21 Mai '13

21 Mai '13

Bom dia a todos, Para recomendar adubação fosfatada é necessário classificar a quantidade de fósforo da análise em conjunto a textura do solo. Tabela da quinta aproximação Nível de Fósforo mg/dm3 Argila% muito baixo Baixo Médio Bom Muito Bom 60-100 <2.7 (menor e igual) 2.8-5.4 5.5-8 8.1-12 >12 35-60 <4 4.1-8 8.1-12 12.1-18 >18 15-35 <6.6 6.7-12 12.1-20 20.1-30 >30 0-15 <10 10.1-20 20.1-30 30.1-45 >45 Fiz o seguinte: solo=c(1:50) # resultado da análise do solo argila=c(60,35,15,0) # teor de argila p1=c(0,2.7,5.4,8,12) # fósforo com + de 60% de argila p2=c(0,4,8,12,18) # fósforo com 35 a 60% de argila p3=c(0,6.6,12,20,30) # fósforo com 15 a 35% de argila p4=c(0,10,20,30,45) # fósforo com < 15% de argila res=c("Muito.Baixo","Baixo", "Medio", "Alto", "Muito.alto") Assim se o resultado de P é 3.5 : 60% de arg, seria classificado como Baixo 35% de arg. seria classificado como muito baixo andei dando uma olhada na solução do post Uso do ifelse A saida desta função seria casada com uma recomendação, Grato a todos -- Hélio Gallo Rocha IFSULDEMINAS - Câmpus Muzambinho

6 8

Ajuda: HPD em conjunto de valores preditos
by Flavia Cristina M.Queiroz 21 Mai '13

21 Mai '13

Olá pessoal, Estou precisando de uma ajuda com o cálculo de intervalo de credibilidade: Tenho um conjunto de dados (valores preditos) com os quais tenho que calcular o intervalo de credibilidade de predição. Estes valores preditos foram calculados em uma implementação de redes neurais. Daí eu calculei o HPD usando a função HPDinterval(), do pacote coda. Fiz assim: despad <- read.table('vintemaio_pad.txt',h=F) dados <- as.matrix(despad) a=density(dados,bw=0.45,from = -1, to = 1) library(coda) arq<-as.mcmc(a$x) mcmcsamp(a$x) HPD<-HPDinterval(arq) HPD Contudo, meu intervalo é muito amplo.. tipo, meus dados vão de -1 a 1 e meu HPD deu [-0.9882583, 0.9138943]. Porque meu intervalo deu tão amplo?? O que teria que fazer de diferente?? ------------------------------------------------------------------------------- Aguardo alguma dica ou solução!! Att., Flávia Cristina

2 2

MOdelo misto sem o grupo para efeito aleatorio
by Pedro Emmanuel Alvarenga Americano do Brasil 21 Mai '13

21 Mai '13

Amigos de R, Estamos tentando implementar em código R uma rotina para suavizar uma curva ROC que utiliza modelos mistos (Haobo Ren , Xiao-Hua Zhou & Hua Liang (2004): A Flexible Method for Estimating the ROC Curve, Journal of Applied Statistics, 31:7, 773-784). O problema é que na hora de rodar a função lme erros diversos aparecem. Aparentemente não estamos sabendo colocar os argumentos de forma apropriada na função, pois o texto original não especifica uma variável de agrupamento. Os objetos são: Rempi <- c(0.75, 0.80, 0.85, 0.85, 0.85, 0.85, 0.85, 0.85, 0.85, 0.85, 0.85, 0.90, 0.95, 0.95, 0.95, 1.00, 1.00, 1.00, 1.00) X <- matrix(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 0.0025, 0.01, 0.0225, 0.04, 0.0625, 0.09, 0.1225, 0.16, 0.2025, 0.25, 0.3025, 0.36, 0.4225, 0.49, 0.5625, 0.64, 0.7225, 0.81, 0.9025, 0.000125, 0.001, 0.003375, 0.008, 0.015625, 0.027, 0.042875, 0.064, 0.091125, 0.125, 0.166375, 0.216, 0.274625, 0.343, 0.421875, 0.512, 0.614125, 0.729, 0.857375),nrow=19,byrow=F) Z <- matrix(c(0, 0, 2e-04, 0.0012, 0.0037, 0.0086, 0.0166, 0.0284, 0.0447, 0.0664, 0.0942, 0.1288, 0.171, 0.2214, 0.281, 0.3504, 0.4304, 0.5217, 0.625, 0, 0, 0, 0, 3e-04, 0.0015, 0.0045, 0.0099, 0.0186, 0.0313, 0.0486, 0.0715, 0.1005, 0.1366, 0.1804, 0.2326, 0.2941, 0.3655, 0.4477, 0, 0, 0, 0, 0, 0, 4e-04, 0.002, 0.0054, 0.0114, 0.0208, 0.0343, 0.0527, 0.0768, 0.1072, 0.1447, 0.1901, 0.2441, 0.3075, 0, 0, 0, 0, 0, 0, 0, 0, 6e-04, 0.0025, 0.0063, 0.013, 0.0231, 0.0376, 0.0571, 0.0823, 0.1141, 0.1531, 0.2002, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1e-04, 9e-04, 0.003, 0.0074, 0.0147, 0.0257, 0.0411, 0.0616, 0.0881, 0.1213, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2e-04, 0.0012, 0.0037, 0.0086, 0.0166, 0.0284, 0.0447, 0.0664, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3e-04, 0.0015, 0.0045, 0.0099, 0.0186, 0.0313, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4e-04, 0.002, 0.0054, 0.0114, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6e-04, 0.0025, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1e-04), nrow=19,byrow=F) grouping <- 1:19 # isto é o que não faz sentido... só escrevi para tentar rodar Em que a expressão do modelo é: Rempi = X*Beta + Z*b + e Tentamos: > lme(fixed = Rempi ~ X, random = ~ Z|grouping) Erro em MEEM(object, conLin, control$niterEM) : Singularity in backsolve at level 0, block 1 Além disso: Mensagens de aviso perdidas: In lme.formula(fixed = Rempi ~ X, random = ~Z | grouping) : Fewer observations than random effects in all level 1 groups Alguma luz pra poder fazer esse modelo dar certo? Abraço forte, Dr. Pedro Emmanuel A. A. do Brasil Curriculum Lattes: http://lattes.cnpq.br/6597654894290806 ResearchGate.net: https://www.researchgate.net/profile/Pedro_Brasil2/ Instituto de Pesquisa Clínica Evandro Chagas Fundação Oswaldo Cruz Rio de Janeiro - Brasil Av. Brasil 4365, CEP 21040-360, Tel 55 21 3865-9648 email: pedro.brasil(a)ipec.fiocruz.br email: emmanuel.brasil(a)gmail.com ---Apoio aos softwares livres www.zotero.org - gerenciamento de referências bibliográficas. www.broffice.org ou www.libreoffice.org - textos, planilhas ou apresentações. www.epidata.dk - entrada de dados. www.r-project.org - análise de dados. www.ubuntu.com - sistema operacional

1 0

Painel com o base
by Alisson Lucrecio 21 Mai '13

21 Mai '13

Caro Colegas, Bom dia. Eu estou precisando montar painel com gráfico feitos no base, vocês poderiam me sugerir a melhor forma de fazer isso? Obrigado. Alisson Lucrécio da Costa

2 1

Acesso ao script da função xyplot do lattice
by Alisson Lucrecio 21 Mai '13

21 Mai '13

Caro Colegas,| Algum de vocês tem acesso ao script da função xyplot do lattice e poderia me envia-lo? Obrigado Alisson Lucrécio da Costa

2 2

ifelse com várias condições
by Fátima Lima Paula 19 Mai '13

19 Mai '13

Prezados, tenho um banco de dados como o abaixo. nome=c("RICARDO CUNHA DA SILVA FILHO","MARIA DE ANDRADE","TEREZA PEREIRA FILHA", "PEDRO ALVARENGA JR","JOSE BRAGANÇA ALMEIDA JR", "GERMANO LIMA BARROS") idade=c(68,67,89,78,65,67) ano=c(1944,1945,1923,1934,1947,1945) dataint=c(20121004,20120907,20120703,20120301,20120805,20121209) banco1=cbind(nome,sobrenome,idade,ano,dataint) banco1=data.frame(banco1) banco1 Separei os nomes em nome1 (primeiro nome), nome2 (nome do meio) e nome3 (último nome), dessa forma: nome=banco1$nome nome=gsub(" D(A|E|O)S? ", " ", nome) separado=strsplit(as.character(nome), split=" ") head(separado) fa=sapply(separado, function(i){ n=length(i) if(n>2) return(c(i[1], paste(i[2:(n-1)],collapse=" "), i[n])) if(n==2) return(c(i[1], NA, i[n])) if(n==1) return(c(i[1], NA, NA)) }) nomenovo=t(fa) head(nomenovo) is.data.frame(nomenovo) nomenovo=as.data.frame(nomenovo) dim(nomenovo) colnames(nomenovo)=c("nome1","nome2","nome3") banco1=cbind(banco1,nomenovo) Depois, vi que existem os apêndices dos nomes (Filho, Filha, Jr etc) Preciso separar os apêndices em uma coluna diferente (nome4). Pensei em criar uma coluna 4 em que tivesse os apêndices do nome3 e no restante viesse NA. Não sei fazer isso. Ainda assim, dessa forma não sei que solução dar para a coluna 3. Pois essa, deveria ter os apêndices substituídos pelos últimos sobrenomes. Minha ideia é: 1)separar o nome2 em 3 de novo assim: nome=banco1$nome2 nome=gsub(" D(A|E|O)S? ", " ", nome) separado=strsplit(as.character(nome), split=" ") head(separado) fa=sapply(separado, function(i){ n=length(i) if(n>2) return(c(i[1], paste(i[2:(n-1)],collapse=" "), i[n])) if(n==2) return(c(i[1], NA, i[n])) if(n==1) return(c(i[1], NA, NA)) }) nomenovo=t(fa) head(nomenovo) is.data.frame(nomenovo) nomenovo=as.data.frame(nomenovo) dim(nomenovo) colnames(nomenovo)=c("nome6","nome7","nome8") banco1=cbind(banco1,nomenovo) 2)Onde for FILHO, JR etc no nome 3, substituir pelo nome8 quando ele não for NA. Pois quando ele é NA, tem que substituir pelo nome6 Não sei se tem como fazer isso. São muitas condições: se for FILHO em nome3 e NA em nome8, substituir por nome6, se for FILHO em nome3 e != de NA em nome8 substituir por nome8. E o pior, tem que fazer isso para FILHO, FILHA, JR etc Alguém sabe fazer isso ou tem outra solução para o problema? Obrigada Fátima -- "Minha felicidade depende da qualidade dos meus pensamentos"

1 0

Indicação de valores fora do intervalo representado em um gráfico.
by Wagner Arbex 17 Mai '13

17 Mai '13

Olá, meus caros, boa tarde; Por favor, gostaria de saber se algum colega poderia me ajudar a resolver o prob a seguir e que não me parece muito difícil, mas não achei uma solução. O exemplo mínimo a seguir mostra cinco conjuntos de dados para desenho de gráficos boxplot no intervalo [-3,3] (essa é a restrição que cria o problema). boxplot(c(1,2,2.6),c(4,5,6),c(-1.5,-3.5,-5),c(1.5,2.5,3.5),c(-2.1,0.2,1),ylim=c(-3,3),horizontal=T) O primeiro e o quinto conjuntos de dados não apresentam nenhum prob, pois seus valores, (1,2,2.6) e (-2.1,0.2,1), estão compreendidos no intervalo [-3,3]. O quarto apresenta seu último valor (3.5) fora do intervalo, mas a representação de todo o quarto conjunto, ( 1.5,2.5,3.5 ), não é problema. O mesmo ocorre com o terceiro conjunto de dados, (1.5,2.5,3.5), dois dos seus valores estão fora do intervalo, mas sua representação não é problema. O problema ocorre qdo TODOS os valores do conjunto estão fora do intervalo [-3,3], como o segundo conjunto, (4,5,6). Nesse caso, o conjunto de dados não é representado na imagem. Meu prob é a necessidade de indicar que esses valores não estão no intervalo que aparece no gráfico. P. ex., para esse caso do segundo conjunto de dados, se aparecesse uma pequena seta, na segunda linha (onde o gráfico seria desenhado), próxima à margem direita do gráfico, apontando para direita, já faria essa indicação. Desde já, agradeço qq sugestão nesse sentido. []s e obrigado. -- Wagner Arbex, DSc Bioinformática e Modelagem Matemática e Computacional de Biossistemas Embrapa Gado de Leite Empresa Brasileira de Pesquisa Agropecuária (Embrapa) Juiz de Fora/MG wagner.arbex(a)embrapa.br Tel.: +55 (32) 3311-7459 | Fax: +55 (32) 3311-7401 www.cnpgl.embrapa.br ____________________________________________________________________________ Aviso de confidencialidade Esta mensagem da Empresa Brasileira de Pesquisa Agropecuaria (Embrapa), empresa publica federal regida pelo disposto na Lei Federal no. 5.851, de 7 de dezembro de 1972, e enviada exclusivamente a seu destinatario e pode conter informacoes confidenciais, protegidas por sigilo profissional. Sua utilizacao desautorizada e ilegal e sujeita o infrator as penas da lei. Se voce a recebeu indevidamente, queira, por gentileza, reenvia-la ao emitente, esclarecendo o equivoco. Confidentiality note This message from Empresa Brasileira de Pesquisa Agropecuaria (Embrapa), a government company established under Brazilian law (5.851/72), is directed exclusively to its addressee and may contain confidential data, protected under professional secrecy rules. Its unauthorized use is illegal and may subject the transgressor to the law's penalties. If you are not the addressee, please send it back, elucidating the failure.

3 2

construção de um objecto para calcular os outliers
by alanarocha＠sapo.pt 17 Mai '13

17 Mai '13

Boa tarde estou a criar objectos para poder fazer o calculo dos outliers mas está a dar tudo NULL, eu queria saber porque.. dados lidos: 2008 2009 2010 2011 2012 1 49 54 57 54 56 2 53 54 55 53 59 3 52 51 52 53 52 4 57 59 59 60 63 5 63 60 63 64 62 6 58 57 59 58 63 estatísticas: #estatísticas descritivas porano<-apply(dados[-c(19,23,36,37),-1],2,function(x)c(summary(x),0)[1:7]) desviopadrao<-apply(dados[-c(19,23,36,37),-1],2,sd,na.rm=TRUE) variancia<-apply(dados[-c(19,23,36,37),-1],2,var,na.rm=TRUE) r<-apply(dados[-c(19,23,36,37),-1],2,range,na.rm=TRUE) diffr<-diff(r) iqr<-apply(dados[-c(19,23,36,37),-1],2,function(x){quantile(x,na.rm=TRUE)[4]-quantile(x,na.rm=TRUE)[2]}) besq<-apply(dados[-c(19,23,36,37),-1],2,function(x){quantile(x, na.rm=TRUE)[2]})-1.5*iqr bdt<-apply(dados[-c(19,23,36,37),-1],2,function(x){quantile(x, na.rm=TRUE)[4]})+1.5*iqr outliers<-list() for(i in names(dados[-c(19,23,36,37),-1])){ tmp<-dados[-c(19,23,36,37),i] outliers[[i]]<-tmp[tmp > bdt[i] | dados[-c(19,23,36,37),i] < besq[i]] } row.names(diffr) <- "Diferen?a do range" # Como o resultado do comando diff ? uma matriz, voc? precisa colocar o nome da linha antes de juntar todos os resultados rownames(porano)<-c("Mínimo","Primeiro Quartil","Mediana","Média","Terceiro Quartil","Máximo","NAs") todos_resultados_dados<-rbind(porano,"Desvio padr?o"=desviopadrao, "Vari?ncia"=variancia, diffr, iqr,besq,bdt,outliers) todos_resultados_dados porano resultado: 2008 2009 2010 2011 2012 Mínimo 48 52 56 55.33 58 Primeiro Quartil 4.708149 4.198439 4.554801 5.255589 4.6888 Mediana 22.16667 17.62689 20.74621 27.62121 21.98485 Média 19 19 19 22 18 Terceiro Quartil 6 5 7 6 7 Máximo 43 44.5 41.5 44 41.5 NAs 67 64.5 69.5 68 69.5 Desvio padr?o Numeric,0 44 Numeric,0 43 Numeric,0 Vari?ncia NULL NULL NULL NULL NULL Diferen?a do range NULL NULL NULL NULL NULL iqr NULL NULL NULL NULL NULL besq NULL NULL NULL NULL NULL bdt NULL NULL NULL NULL NULL outliers NULL NULL NULL NULL NULL porquê NULL? Cumprimentos Ana Rocha

2 1