Maio 2012 - R-br - DINF UFPR

Re: [R-br] Teste de médias em GLMM com covariável!
by Ivan Bezerra Allaman 10 Mai '12

10 Mai '12

Valeu Walmes, era exatamente isso que eu procurava. Agora entendi. Muito obrigado! \begin{signature} <<>>= Prof. Dr. Ivan Bezerra Allaman Universidade Estadual de Santa Cruz Departamento de Ciências Exatas e Tecnológicas Ilhéus/BA - Brasil Fone: +55 73 3680-5596 E-mail: ivanalaman@yahoo.com.br/ivanalaman@gmail.com @ \end{signature} ________________________________ De: Walmes Zeviani <walmeszeviani(a)gmail.com> Para: Ivan Bezerra Allaman <ivanalaman(a)yahoo.com.br> Enviadas: Quinta-feira, 10 de Maio de 2012 0:14 Assunto: Re: [R-br] Teste de médias em GLMM com covariável! Ivan, Bem, de fato não estou seguro sobre ter entendido o que você quer, mas vamos lá. Partindo do fato que o nosso modelo aponta interação entre F e x nos teremos retas não paralelas (na escala do preditor linear), portando diferenças de médias não são constantes a medida que variamos o valor de x. Podemos comparar as médias em um valor fixo de x. Por exemplo, podemos comparar o peso final de suínos do sexo macho vs fêmea (F) supondo que os animais tem mesma idade (x=100 dias). Outra possibilidade, que deve ser a que você quer, é comparar nas médias amostrais de x nos níveis de F. Assim, comparamos o peso final dos machos com idade média de 100 dias com as fêmeas com idade média final de 95 dias. Se isso é uma comparação justa ou não é outra história e depende do contexto. Eu fiz o CMR desse último caso. É isso que você procura? Eu prefiro sempre que possível evitar as comparações (principalmente seguidas de letras) e dar um gráfico com curvas ajustadas acompanhadas das bandas de confiança. Assim o leitor pode ver como são as diferenças em cada nível de x, entender a relação de x com F e apontar para quais intervalos em x existe ou se pronuncia a diferença entre os níveis de F. É sem dúvida mais informativo. Espero ter ajudado. pred <- expand.grid(F=levels(da$F), x=seq(0,1,l=30)) pred$eta <- predict(ac0, newdata=pred) pred$y <- predict(ac0, newdata=pred, type="response") require(lattice); require(latticeExtra) # retas não paralelas. em que ponto vou compara-las? # em cada lugar dá uma conclusão diferente! xyplot(eta~x, groups=F, data=pred, type="l") # você pode comparar todas em um valor fixo, como x=0.5 # ou pode comparar ela no valor médio de x para cada nível de F xm <- with(da, tapply(x, F, mean)); xm mpm==0.5 mpm2 <- mpm for(i in 1:length(xm)) mpm2[i,mpm2[i,]==0.5] <- xm[i] mpm2 mpm2%*%betas exp(mpm2%*%betas) # diferença nas médias populacionais marginais duas à duas nlevels(da$F) comp <- t(combn(1:nlevels(da$F), 2)) rownames(comp) <- apply(comp, 1, paste, collapse="-") dmpm <- t(apply(comp, 1, function(vs) mpm2[vs[1],]-mpm2[vs[2],])) rownames(dmpm) <- rownames(comp) apply(dmpm, 1, function(x) x%*%betas) require(multcomp) summary(glht(ac0, linfct=mpm2)) # médias populacionais marginais summary(glht(ac0, linfct=dmpm)) # contraste entre as médias populacionais marginais xyplot(y~x, groups=F, data=da, type="p")+ as.layer(xyplot(y~x, groups=F, data=pred, type="l")) À disposição. Walmes. ========================================================================== Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W) Departamento de Estatística - Universidade Federal do Paraná fone: (+55) 41 3361 3573 VoIP: (3361 3600) 1053 1173 e-mail: walmes(a)ufpr.br twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================

1 0

Teste de médias em GLMM com covariável!
by Ivan Bezerra Allaman 10 Mai '12

10 Mai '12

Boa tarde senhores! Estou a alguns dias pesquisando em como fazer um teste de comparações em modelo com covariáveis e cujos slopes são diferentes significativamente, já que, as comparações devem ser feitas com as médias ajustadas para as covariáveis. Já procurei no nablle e pela internet, mais nada esclarecedor. Se alguém já trabalhou com este tipo de problema e que queira me dar uma dica, principalmente no que diz respeito a algum material que explique filosoficamente e matricialmente como funciona as comparações múltiplas neste casos de covariável e ficarei extremamente grato. Segue um CRM do que eu vos falo. set.seed(1) da <- data.frame(Trat=gl(5,10),covi=sample(rpois(100,30),50,replace=TRUE),vari=sample(rpois(100,60),50,replace=TRUE),animal=seq(1:50)) da library(lme4) ancova <- glmer(vari ~ Trat*covi + (1|animal%in%Trat),data=da,family=poisson) Anova(ancova)#A interação mostra que as diferenças de slopes entre Trat são significativas ancova2 <- update(ancova,~.-Trat:covi) anova(ancova2,ancova) library(multcomp) summary(glht(ancova,linfct=mcp(Trat='Tukey'),interaction_average=TRUE,covariate_average=TRUE))#Acredito que isto é feito sem correção para as covariáveis. Em algumas buscas, encontrei que a função effect da library ``effects'' faz tal ajutes, mais não estou seguro quanto ao seu uso. Por isso, com o entendimento de como funciona isto matricialmente e acho que consigo resolver via 'glht'. \begin{signature} <<>>= Prof. Dr. Ivan Bezerra Allaman Universidade Estadual de Santa Cruz Departamento de Ciências Exatas e Tecnológicas Ilhéus/BA - Brasil Fone: +55 73 3680-5596 E-mail: ivanalaman@yahoo.com.br/ivanalaman@gmail.com @ \end{signature}

4 5

Bandas de confiança em plot de dispersão
by Marcelo Cardoso mello 09 Mai '12

09 Mai '12

Pessoal, gostaria de adicionar bandas de confiança na CMR que segue. Agradeço qualquer dica. x<-structure(list(week = c(9L, 10L, 11L, 12L, 9L, 10L, 11L, 12L, 9L, 10L, 11L, 12L, 9L, 10L, 11L, 12L, 9L, 10L, 11L, 12L, 9L, 10L, 11L, 12L, 9L, 10L, 11L, 12L, 9L, 10L, 11L, 12L, 9L, 10L, 11L, 12L, 9L, 10L, 11L, 12L, 9L, 10L, 11L, 12L, 9L, 10L, 11L, 12L, 9L, 10L, 11L, 12L, 9L, 10L, 11L, 12L, 9L, 10L, 11L, 12L, 9L, 10L, 11L, 12L, 9L, 10L, 11L, 12L), kid = c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 14L, 14L, 14L, 14L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 17L, 17L, 17L, 17L), treatment = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L ), gain = c(1.2, 1, 1.1, 1.3, 1.2, 1.1, 1.4, 1.5, 1.3, 1.4, 1.4, 1.6, 1.1, 1.1, 1.2, 1.3, 1.2, 1.3, 1.2, 1.3, 1.1, 1.1, 1.1, 1.2, 1.1, 1.2, 1.3, 1.5, 1.3, 1.3, 1.3, 1.4, 1.2, 1.5, 1.9, 2.1, 1.3, 1.2, 1.4, 1.7, 1.5, 1.7, 1.6, 1.7, 1.4, 1.5, 1.7, 1.8, 1.2, 1.2, 1.4, 1.6, 1, 1.1, 1.4, 1.5, 1.4, 1.8, 2.1, 2.1, 1.1, 1.3, 1.4, 1.8, 1.2, 1.5, 1.7, 1.9)), .Names = c("week", "kid", "treatment", "gain"), class = "data.frame", row.names = c(NA, -68L)) y<-aggregate(.~week+treatment, x, FUN=mean) s1<-subset(y, treatment==1) s2<-subset(y, treatment==2) plot(gain~week, data=s1, col=1, type="o", ylim=c(1,2), lty=1) points(gain~week, data=s2, col=2, type="o", lty=2) legend(9.5,1.7, c("t1","t2"), col=c(1,2), lty=c(1,2))

3 2

Dúvida urgente sobre array
by Victor Eduardo 09 Mai '12

09 Mai '12

Galera, queria tirar uma dúvida sobre os arrays, isso está me dando dor de cabeça para realizar uma função para gerar um modelo. Um exemplo seria: > ar_K<-array(1,dim=c(1,2,2)) > ar_K , , 1 [,1] [,2] [1,] 1 1 , , 2 [,1] [,2] [1,] 1 1 > ar_K[,,1] [1] 1 1 > is.matrix(ar_K[,,1]) [1] FALSE > is.vector(ar_K[,,1]) [1] TRUE Queria selecionar uma matriz desse array, mas sempre vem como vetor. Como posso selecionar só a matriz dentro do array? Abraços!!!

2 3

Label=T em barplot
by Diego Leal 09 Mai '12

09 Mai '12

Boa noite pessoal, estou usando a funcao barplot() para gerar um grafico. Existe algum parametro equivalente ao label=T da funcao hist() que eu possa usar e que mostre os valores plotados? esp=c(0,12,30,5,2,1,0) barplot(rev(esp),horiz=TRUE,names.arg=(length(esp)-1):0) gostaria que os valores 0, 12, 30 ... aparececem rotulando as barras sem ter que faze isso manualmente. Att, Diego Leal Togni

3 2

Duvida sobre a função query() do pacote seqinr
by Augusto Ribas 09 Mai '12

09 Mai '12

Ola pessoal, eu estou tentando procurar e extrair algumas sequências de DNA pelo R Em um livro que estou lendo fui guiado até o pacote seqinr, ele tem uma função chamada query() pra fazer isso. A principio tudo parece bem simples e intuitivo. Por exemplo, eu tenho interesse na sequência 16s rRNA de varias especies de sapinhos. Ai se procuro pra uma especie tenho: > library(seqinr) > s <- choosebank("genbank") > query("SP","SP=Rhinella rubescens AND K=16S RIBOSOMAL RNA",s$socket) > SP$nelem [1] 1 > getName(SP) [1] "GU907196.RR1" Blz eu achei 1 sequência especifica pro que eu quero, e o nomezinho pra acessar ela. Procuro outra > query("SP","SP=Leptodactylus syphax AND K=16S RIBOSOMAL RNA",s$socket) > print(SP$nelem) [1] 2 > print(getName(SP)) [1] "JF789923" "JF789924" Blz, 2 sequencias,ai eu fui ao site do genbank, que é a base da dados da uma olhada e estao la blz So que muitas especies que procuro o R para. Por exemplo: > query("SP","SP=Leptodactylus marmoratus AND K=16S RIBOSOMAL RNA",s$socket) O R fica parado e tenho que fechar ele. Eu não entendo o que esta dando errado, principalmente pq algumas especies funcionam e outras não. Essa especie em especifico num era pra achar nada, não achei registro dela procurando pelo site do genbank. Mas a hora que vou usar o R num vai. Outro exemplo é "Leptodactylus syphax" que devia acha 2 sequencias mas ele para também. E a maior parte das especies num vai. Eu uso windows 7 como SO. Estou com o R 2.15.0, internet residencial. Alguem tem experiencia com isso? Sabe dizer o que esto fazendo de errado? E alguem conhece alguma outra forma de fazer o mesmo tipo de busca usando o R? O task view indica o Biocondutor que tem o pacote SRAdb que deve fazer a mesma coisa, mas não consegui usar. Eu vi que o livro que é de 2011 a função mudou a ordem dos argumentos, e teve um update no pacote, mas não sei se estou fazendo algo errado, se tenho que fazer mais alguma coisa. Os manuais que vi não tem um faq com problemas comum e soluções. Minha lista de especies caso alguem queria ver é: lista.spp<-c("Rhinella rubescens", "Rhinella fernandezae", "Elosia rustica", "Crossodactylus gaudichaudii", "Leptodactylus pustulatus", "Leptodactylus syphax", "Dendropsophus cachimbo", "Leptodactylus marmoratus", "Physalaemus soaresi", "Hylodes nasus", "Scinax fuscomarginatus", "Leptodactylus petersii", "Chiasmocleis capixaba", "Ceratophrys aurita", "Pleurodema diplolister", "Cystignatus gigas", "Scinax trilineatus", "Elosis nasus", "Dryadophis bifossatus" ) Se alguem puder dar uma luz, não sei nem por onde começar a procurar pra entender o que ta acontecendo. Vou quebrando a cabeça por enquanto, Boa noite a todos. -- Grato Augusto C. A. Ribas Site Pessoal: http://augustoribas.heliohost.org Lattes: http://lattes.cnpq.br/7355685961127056

2 4

Criar Buffer
by Ricardo Dobrovolski 09 Mai '12

09 Mai '12

Olá, Como posso criar um buffer a partir de uma linha? Essa linha foi criada a partir de um par de pontos. Eu preciso desse buffer no formato raster. Abraço, RD.

1 0

Como ler formato dbc no R
by Bernardo Rangel Tura 09 Mai '12

09 Mai '12

Primeiro de tudo bom dia a todos Na minha atividade costumo a usar bancos de dados do Ministério da Saúde. Recentemente precisei trabalhar com todas as internações do país disponíveis para análise. O que dava uma quantidade boa de arquivos pois são 17 anos de arquivos mensais por estado = 5.486 arquivos. Mas isto não chega a ser um problema. O problema é que o arquivo vem num formato de dados .dbc, que até onde consegui descobrir é de DBF compactado utilizado pelo foxpro. Como eu uso linux tenho que emular um windows (wine) utilizar o programa tabwin para transformar dbc em dbf. Como isto o conjunto de dados sai de 10 Gb para 35 Gb. No intuito de simplificar minha vida gostaria de saber de vocês; é possível um mortal que utiliza R 2.15 - atualizado e com praticamente todo e qualquer pacote instalado - rodando Ubuntu 12.04 64 bits num phenom 2 X6 ler diretamente estes arquivos .dbc ? Se for possível o que devo fazer? Antes de mais nada envio o menor arquivo que trabalho como Minimo Exemplo Reprodutível. Ele esta em 2 formatos .dbc e .dbf corresponde ao extrato das internações de Roraima de dezembro de 1997 Certo de sua atenção e genialidade já me despeço agradecendo []s Tura

1 0

Tukey para parcelas subdividad
by Gilson Sanchez 08 Mai '12

08 Mai '12

Pessoal, alguém tem um script para análise de dados em DBC com parcelas subdividas, alias preciso do teste de Tukey (5% de prob). Obrigado -- Bach. Agr.Gilson Sánchez Chia, MSc. Laboratório de Fisiologia Vegetal Embrapa Amazônia Ocidental Fone: (92) 3303-7841 / 3303-7919 http://lattes.cnpq.br/3684553005529634 --- Antes de imprimir, pense na sua responsabilidade com o Meio Ambiente. Please consider the environment before printing this email. Renew, Reduce and Recycle. The Planet Knows.*** Salva un Árbol.* * *** * PIENSA**** EN EL MUNDO PARA TUS HIJOS** *

5 6

Contraste
by Adriano Carvalho Costa 08 Mai '12

08 Mai '12

Boa noite, Gostaria de saber se alguém tem um scrip para contraste. Obrigado!

3 2