Julho 2014 - R-br - DINF UFPR

Teste para comparar distribuições
by Luiz Roberto Martins Pinto 24 Jul '14

24 Jul '14

Caros amigos, Fiz uma pesquisa no Google, mas não consegui uma solução satisfatória. Então peço ajuda. Preciso de encontrar um teste que *quantifique* a semelhança entre dois vetores para checar se eles têm o mesmo tipo de distribuição, por meio de p.value. O interesse principal é comparar os vetores x e ab. Data set: http://www.datafilehost.com/d/f3757310 library(Hmisc) library(grDevices) load('NNC D_d3kk2.RData') # Avaliação qualitativa (gráfica): dev.set(which=1) bpplot(d3) boxplot(d3) # Aparentemente os vetores x, abC e ab são muitos semelhantes # A minha expectativa é encontrar um teste que quantifique esta semelhança por meio de p.value # Então fiz os testes abaixo # Avaliação quantitativa: v2=d3$ab # Comparando o vetor x com o vetor ab: v1=d3$x round(ks.test(v1, v2)$p.value,2) # Out put: p.value=0. Então o vetor x é estatisticamente DIFERENTE do vetor ab. # Comparando o vetor abC com o vetor ab: v1=d3$abC round(ks.test(v1, v2)$p.value,2) # Out put: p.value=0.65. Então o vetor x é estatisticamente IGUAL do vetor abC. # Comparando o vetor abC com o vetor x: v2=d3$x round(ks.test(v1, v2)$p.value,2) # Out put: p.value=0. Então o vetor abC é estatisticamente DIFERENTE do vetor x. ### *Testes alternativos*, para comparar o vetor x com ab: v2=d3$ab v1=d3$x t1=round(ks.test(v1, v2)$p.value,2);t1 # Out put: p.value=0 t2=round(ks.test(v1, v2, alternative = "l")$p.value,2);t2 # Out put: p.value=0 t3=round(ks.test(v1, v2, alternative = "g")$p.value,2);t3 # Out put: p.value=0.85 # Então, para este Data set (NNC D_d3kk2.RData), e para o t3=0.85. Então, por meio deste teste, o vetor x é estatisticamente IGUAL do vetor ab. # *No entanto*, fazendo os testes com conjuntos de dados diferentes (10000 Data set) (mas, sempre com mesmo tamanho de vetor), quando comparo o vetor x com o vetor ab: ## em aproximadamente 60% das comparações (ou em 6000 Data set) t1>0.05 ## em aproximadamente 80% das comparações (ou em 8000 Data set) t2>0.05 ## em aproximadamente 80% das comparações (ou em 8000 Data set) t3>0.05 ## em aproximadamente 85% das comparações (ou em 8500 Data set) t1 *OU* t2 >0.05 ## em aproximadamente 85% das comparações (ou em 8500 Data set) t1 *OU* t3 >0.05 ## em 0% das comparações (ou em 0 Data set) t2 *E* t3 >0.05 ## em 100% das comparações (ou em 10000 Data set) t2 *OU *t3 >0.05 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

3 4

Re: [R-br] Como fazer mapa temático
by Waldir De Carvalho Junior 22 Jul '14

22 Jul '14

Prezada Jacqueline sua questão está muito genérica. É preciso ser mais específico para te ajudar. No entanto, os pacotes 'sp', 'rgdal', 'plotkml', entre outros, tem a possibilidade de plotar os mapas. verifique Att __________ Dr. Waldir de Carvalho Junior Engº Agroº D.Sc. Ciência do Solo - Mapeamento Digital Embrapa Solos Empresa Brasileira de Pesquisa Agropecuária (Embrapa) Rio de Janeiro/RJ Waldir.Carvalho(a)embrapa.br Telefone: +55 (21) 2179-4614 +55 (21) 2179-4614 | Fax: +55 (21) 2274-5291 www.embrapa.br/solos | twitter.com/embrapa Confira também: www.facebook.com/agrosustentavel Call Send SMS Add to Skype You'll need Skype Credit Free via Skype ----- Mensagem original ----- De: r-br-request(a)listas.c3sl.ufpr.br Para: r-br(a)listas.c3sl.ufpr.br Enviadas: Terça-feira, 22 de julho de 2014 12:00:02 Assunto: Digest R-br, volume 43, assunto 22 Enviar submissões para a lista de discussão R-br para r-br(a)listas.c3sl.ufpr.br Para se cadastrar ou descadastrar via WWW, visite o endereço https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br ou, via email, envie uma mensagem com a palavra 'help' no assunto ou corpo da mensagem para r-br-request(a)listas.c3sl.ufpr.br Você poderá entrar em contato com a pessoa que gerencia a lista pelo endereço r-br-owner(a)listas.c3sl.ufpr.br Quando responder, por favor edite sua linha Assunto assim ela será mais específica que "Re: Contents of R-br digest..." Tópicos de Hoje: 1. Como fazer mapa temático (Jacqueline Cantú) ---------------------------------------------------------------------- Message: 1 Date: Tue, 22 Jul 2014 17:16:55 +0300 From: Jacqueline Cantú <jacqueline.cantu(a)hotmail.com> To: "r-br(a)listas.c3sl.ufpr.br" <r-br(a)listas.c3sl.ufpr.br> Subject: [R-br] Como fazer mapa temático Message-ID: <BLU177-W17CD82038F1F8B202CBD5A8BF10(a)phx.gbl> Content-Type: text/plain; charset="iso-8859-1" Olá pessoal,Construí o semivariograma cruzado, ajustei modelo e escolhi o melhor modelo e preciso fazer o mapa temático para cada uma das variáveis. Alguém sabe me dizer como fazer o mapa temático? Muio obrigada pela atenção.Att, Jacqueline -------------- Próxima Parte ---------- Um anexo em HTML foi limpo... URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20140722/ad115a93/atta…> ------------------------------ Subject: Legenda do Digest _______________________________________________ R-br mailing list R-br(a)listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br ------------------------------ Fim da Digest R-br, volume 43, assunto 22 ***************************************** ____________________________________________________________________________ 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.

1 0

Identificar pontos num gráfico e adicionar rótulos
by Helder Gramacho 22 Jul '14

22 Jul '14

Bom dia pessoal, Estou tentando usar o identify() para identificar pontos num gráfico e em seguida adicionar um rótulo para os pontos identificados, este assunto já havia sido discutido aqui na lista e foi resolvido com o comando abaixo, no meu caso estou utilizando o RStudio, só que depois que executo o comando parece que o gráfico trava e os rótulos não são adicionados, está faltando algo nos comandos? ### <code r> dados <- data.frame(x=1:10,y=1:10,Nome=paste('ponto',letters[1:10])) plot(dados$x,dados$y) identify(dados$x,dados$y, labels=dados$Nome) ### </code> Obrigado, *Hélder Gramacho * Recife-PE / *agrohelder(a)gmail.com <agrohelder(a)hotmail.com>*

3 4

Selecionar vários níveis dentro de um fator com for()
by ASANTOS 16 Jul '14

16 Jul '14

Boa noite pessoal, Tenho um fator (indiv) no objeto dados_cal[,3], composto de vários níveis sendo soldado1, .... ate soldado 20 e operario1, ... ate operario20, gostaria de selecionar só os operário, para isto eu fiz: #Operários x<- c("operario1","operario2","operario3","operario4","operario5","operario6", "operario7" ,"operario8","operario9","operario10","operario11","operario12", "operario13","operario 14","operario15","operario16","operario17","operario18", "operario19","operario20") RES=NULL for(l in 1:length(x)){# Rayon des cercles teste<-dados_cal[dados_cal[,3]==x[l]] RES=rbind(RES,teste) } # Mas o meu for() aplicado ao objeto x não funcionou, alguém poderia dar um help? 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 2

Mensagens de erro no R
by Clécio da Silva Ferreira 16 Jul '14

16 Jul '14

Meus caros, aconteceu comigo hoje (1ª vez) ao tentar instalar e utilizar o help de funções específicas de um pacote (instalei hoje gamlss). O pacote foi corretamente instalado. Ao solicitar carregar pacotes, mesmo instalados anteriormente, começaram a aparecer mensagens como: library(MASS) Mensagens de aviso perdidas: In fetch(key) : internal error -3 in R_decompress1 ?memory.size Mensagens de aviso perdidas: In fetch(key) : internal error -3 in R_decompress1 Aumentei então o limite de memória: memory.size() [1] 3232.17 > memory.limit() [1] 3975 > memory.size(max=T) [1] 3890.69 > memory.limit(size=10000) [1] 10000 Porém, o help de funções do gamlss não está funcionando corretamente. Ao solicitar o help de alguma função, aparece: Error in fetch(key) : lazy-load database '' is corrupt Atualizei os pacotes, mas não resolveu o problema. Alguma ideia? Atenciosamente, Prof. Clécio da Silva Ferreira Depto de Estatística - UFJF Tel: (32) 2102-3306 Página: www.ufjf.br/clecio_ferreira

4 4

Novas versões dos pacotes R
by Luiz Roberto Martins Pinto 15 Jul '14

15 Jul '14

Caros amigos da lista. Eu sou um iniciante no uso do R. Depois que comecei a utilizá-lo a análise de dados ficou mais divertida e criativa. Sinto-me livre... incrível!!! O uso do Tinn-R também é um fator motivacional para mim, devido à simplicidade de uso. Todavia, como iniciante, tenho tido algumas dificuldades em referência às atualizações do R e dos pacotes. Certa vez atualizei o R e algumas funções/pacotes deixaram de funcionar. Quase entrei em pânico, posto que jamais pensei que a simples atualização do R pudesse causar tanto transtorno. Acabei gastando muito tempo procurando erros no script, até que alguém sugeriu que o script poderia não rodar pela simples atualização do R. A atualização de pacotes é outro problema. Quando atualizei o R na última vez uma função que utilizava parou de funcionar. Então precisei de refazer parte importante do script. Diante do exposto vejo que preciso ter: 1) cópia dos pacotes; 2) cópia da versão do R utilizada em cada script É assim que vocês procedem? É possível que o R seja melhorado/atualizado sem prejuízo no uso dos scripts? É possível que os autores dos pacotes façam suas atualizações sem perda de funcionalidade das funções anteriormente liberadas? Ex. liberar o LSmatrix e LSmeans (informando que esta [e uma versão mais atualizada) e manter o popMeans e popMatrix? É querer muito? Talvez, para o expert em R estas questões sejam simplórias, e até idiotas, mas para alguns poucos (espero) simples usuários como eu estas questões são relevantes. Luiz Roberto. 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

5 10

modelagem (com glmer) depois de scale() de variáveis contínuas
by Pedro Emmanuel Alvarenga Americano do Brasil 12 Jul '14

12 Jul '14

Amigos de R, Eu estou tentando fazer uma modelagam logistica multinível com com alguns modelos especificados a priori. São sete ao total. Em 5 deles aparece o seguinte aviso. > ptm <- Sys.time() > f6a <- glmer(Hosp.Death ~ Adm.Diag.Class + SAPS3PT + Chron.Health.Status.Rev + LogNDiasantesICU + V2.Total.interact + Accredited.YN + Resid.UTI + MedSurg.ICU + ICUbed1020 + HCP.Cover.Class + (1 | ICUCode), data = b, family = binomial,nAGQ=1) # Warning message: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue - Rescale variables? > Sys.time() - ptm Time difference of 14.5636 mins Infelizmente o modelo não diz quais as variáveis que estão atrapalhando mas eu suspeito de duas das quatro que há no modelo. Eu coloquei esse erro no google e apareceram alguns comentários do desenvolvedor no StackOverflow que eu achei interessantes. Na essência, ele faz uma recodificação sugerida pelo aviso com a função base::scale() que centraliza a distibuição das variáveis continuas dessa forma: > z <- b > z[,numcols] <- scale(z[,numcols]) > summary(z[,numcols]) Age SOFAPT SAPS3PT HospDayPriorICU Min. :-2.4008 Min. :-0.7762 Min. :-2.5456 Min. : -0.06196 1st Qu.:-0.6911 1st Qu.:-0.7762 1st Qu.:-0.6716 1st Qu.: -0.06196 Median : 0.1378 Median :-0.4487 Median :-0.1362 Median : -0.06196 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000 Mean : 0.00000 3rd Qu.: 0.8113 3rd Qu.: 0.5336 3rd Qu.: 0.6000 3rd Qu.: -0.03865 Max. : 2.5210 Max. : 6.0999 Max. : 5.9542 Max. :222.87390 > summary(b[,numcols]) Age SOFAPT SAPS3PT HospDayPriorICU Min. : 16.00 Min. : 0.00 Min. : 5.00 Min. : 0.000 1st Qu.: 49.00 1st Qu.: 0.00 1st Qu.: 33.00 1st Qu.: 0.000 Median : 65.00 Median : 1.00 Median : 41.00 Median : 0.000 Mean : 62.34 Mean : 2.37 Mean : 43.03 Mean : 2.653 3rd Qu.: 78.00 3rd Qu.: 4.00 3rd Qu.: 52.00 3rd Qu.: 1.000 Max. :111.00 Max. :21.00 Max. :132.00 Max. :9564.000 Em seguida fazer um update do modelo com os novos dados. O que remove os avisos e > f6az <- update(f6a,data=z) Olhando o summary dos dois modelos, na essência é o mesmo modelo com o mesmo AIC, BIC, com a mesma variação do efeito aleatório e as mesma variáveis significativas no efeito fixo. Na verdade então eu nem sei muito bem porque o autor recomendaria essa transformação. Mas agora que eu quero fazer umas análises de previsões medias por grupos marginalizando a variação do efeito aleatório, a pergunta é como fazer a destranformação do scale(), ou seja, qual seria a correspondência dos valores dados pelo scale e os valores originais da escala da variável? Dei uma olhada na documentação da função scale() e entendi que essa centralização ocorre subtraindo-se os valores da média e depois dividindo-se pelo sd(). Assim, entendo que a "destroanformação" seria multiplicando-se pelo sd() e adicionando-se a média. Dessa forma bateu com os valores originais. > as.vector(scale(b$SAPS3PT[1:30])) * sd(b$SAPS3PT[1:30]) + mean(b$SAPS3PT[1:30]) [1] 64 37 59 42 32 51 69 58 55 28 29 65 43 54 51 50 52 57 16 56 46 56 34 31 81 [26] 85 61 80 26 80 > b$SAPS3PT[1:30] [1] 64 37 59 42 32 51 69 58 55 28 29 65 43 54 51 50 52 57 16 56 46 56 34 31 81 [26] 85 61 80 26 80 Como eu nunca fiz isso antes, eu queria saber se voces acham que essa centralização das variáveis contínuas é razoável, se a destransformação é assim mesmo e se fariam essa análise de forma diferente? Abraços, Dr. Pedro Emmanuel A. A. do Brasil http://blog.ipec.fiocruz.br/lapclin-chagas/ Curriculum Lattes: http://lattes.cnpq.br/6597654894290806 ResearchGate.net: https://www.researchgate.net/profile/Pedro_Brasil2/ Instituto Nacional de Infectologia/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 e-mail: pedro.brasil(a)ipec.fiocruz.br e-mail: emmanuel.brasil(a)gmail.com

2 1

Caret
by Gleyson Vieira dos Santos 11 Jul '14

11 Jul '14

Nome: Gleyson Gostaria de saber como consigo calcular o quadrado médio do resíduo (QMR), e como obtenho os valore preditos? -- GLEYSON VIEIRA DOS SANTOS ZOOTECNISTA - UFPI DOUTORANDO EM CIÊNCIA ANIMAL-UFPI TEL: (89) 9985-5488

2 1

Recomendação de adubação
by Hélio Gallo Rocha 11 Jul '14

11 Jul '14

Caros listeiros e interessados em recomendação de adubação. Estou precisando de ajuda para realizar recomendação de adubação. Como exemplo utilizei a cultura do café com cálculos realizados de acordo com a quinta aproximação para Minas Gerais, que pode ser obtido no seguinte link, que no caso, é especifico para a cultura do café, e a tabela em questão se encontra na pag. 4. http://www.dpv24.iciag.ufu.br/new/dpv24/Apostilas/Cafe%20-%205%20Aproximaca… Uma parte já está pronta, postada em msg anterior pelo mestre Elias. 1. Classifica-se o teor de fósforo em cinco classes, de acordo com o teror (%) de argila, segundo este CRM: ## Dados p <- c( 7, 8 , 7, 7, 10);p arg <- c(10 , 20, 40, 70, 70); arg ## Classificação da argila class_arg=c(0,15.01,35.01,60.01,100) nivel_arg=c("Arenoso","Siltoso","Franco_Argiloso","Argiloso") clas_arg=cut(arg,class_arg,nivel_arg) ## Visualizando a classificação das amostras de argila data.frame(arg,clas_arg) # Classificação do Fósforo de acordo com nível de argila escada <- function(x, y, x.br, y.br, ...) { ix <- findInterval(x, x.br, ...)+1 lr <- lapply(1:nrow(y.br), function(j) { i <- which(ix==j) list(i=i, r=findInterval(y[i], y.br[j,], ...)+1) }) i <- sort(unlist(lapply(lr, function(x) x$i)), index.return=TRUE, method='quick')$ix unlist(lapply(lr, function(x) x$r))[i] } ## teores de fosforo (mg/dm3) para a cultura do café que serão classificados ## de acordo com quatro classes de solo de acordo com teror (%) de argila cafefosf <- rbind(c(7.5,15.1,22.6,33.81),c(5.01,9.1,15.1,22.51), c(3.1,6.1,9.1,13.51 ), c(2,4.1,6.1,9.01)) ## Determinando a classe de ferfilidade(cf.cafe) do solo para a cultura do café cf.cafe=escada(arg, p, c(15.001, 35.001, 60.001),cafefosf) ## verificando a classe de fertilidade de acordo com quantidade de fósforo ##(mg/dm3) e teor de argila no solo, onde arg = %arg, p = fósforo, cf.cafe = classe ## é um exemplo com cinco classes cbind(arg,p,cf.cafe) 2. Para fazer a recomendação de adubação é realizada de acordo com a produção esperada(sacas café benificiado ha-1), segundo *produ.cafe* produ.cafe <- c(0,20,30,40,50,60,100) 3. O resultado é de acordo com *aduba.cafe* ( em kg/ha/ano) aduba.cafe <- rbind(c(80,60,40,20,0),c(70,55,35,18,0),c(60,50,30,15,0),c(50,40,25,0,0),c(30,20,10,0,0));aduba.cafe Para verificar se deu certo, a primeira linha do resultado de cbind(cf.cafe) que é 1, com uma produção estimada em 25 sacas daria como resultado a recomendação de 30 kg/ha/ano de fósforo. Agradeço desde já. -- Hélio Gallo Rocha IFSULDEMINAS - Câmpus Muzambinho

1 0

popMatrix
by ana paula coelho madeira 10 Jul '14

10 Jul '14

Boa tarde! Prezados colegas, a função popMatrix foi removida do pacote doBy. Como posso substituí-la? Agradeço, Ana Paula

2 1