- R-br - DINF UFPR

Ler dados do epi info
by Emerson Cotta Bodevan 11 Mai '12

11 Mai '12

Prezados, como faço para importar dados do epi info para o R.? -- *Emerson*

4 4

Simulação de um DQL no R
by Edmundo Caetano 11 Mai '12

11 Mai '12

Prezados, Pretendo simular um experimento em quadrados latinos no R. Como devo proceder??? Edmundo Caetano

4 8

Problemas com comando witch
by Tito Conte 11 Mai '12

11 Mai '12

tenho dois vetores: p=c(2,3,4,5,6,7,8,9) r=c(3,4,5) quando uso which(p==r) ou which(r==p) ele me retorna as posições onde r ocorre em p como fazer para saber onde p ocorre em r? obrigado (já dei um google nisso e não encontrei ) Tito Conte

4 7

Escalas
by Eric Keven 10 Mai '12

10 Mai '12

Boa noite, Como eu faço para alterar a escala dos gráficos e fazer com que o programacomece num valor que eu fixar para a escala?Obrigado! Att, Eric Keven Silva Coordenador do Núcleo de Estudos de Processos da Indústria de Alimentos (PROALI) Graduando em Engenharia de Alimentos Universidade Federal de Lavras

5 4

Comparações múltiplas usando MMC (mean–mean multiple comparisons) plots
by ASANTOS 10 Mai '12

10 Mai '12

*Bom dia pessoal, Estou tendo problemas em visualizar os resultados de uma comparação múltipla usando GLM, com o uso MMC plots do pacote HH. Bom, tenho quatro grupos funcionais de cupins que são humívoros (humi), ceifadores(ceifa), intermediários (inter) e xilófagos (xilo) que são meus tratamentos e como variável resposta o número de ocorrência de espécies de cada grupo (acorr). Para coleta da variável resposta, fizemos 6 transectos que são as repetições. Para ver se existe diferença entre os grupos funcionais fiz:* # mod2<-glm(ocorr~grupo, family="quasipoisson", data=trat1_90antes)## Um plot(mod2) prévio mostrou que quasipoisson foi a distribuição mais adequada aos dados que poisson summary(mod2) glm(formula = ocorr ~ grupo, family = "quasipoisson", data = trat1_90antes) Deviance Residuals: Min 1Q Median 3Q Max -1.9148 -0.4946 0.1213 0.4178 1.3813 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.5041 0.1624 9.260 1.13e-08 *** grupohumi -1.0986 0.3249 -3.382 0.00296 ** grupointer -0.8979 0.3019 -2.974 0.00750 ** grupoxilo -0.6568 0.2780 -2.363 0.02838 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for quasipoisson family taken to be 0.7123377) Null deviance: 28.663 on 23 degrees of freedom Residual deviance: 16.889 on 20 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 5 anova(mod2, test="Chi") ## Anova do modelo Analysis of Deviance Table Model: quasipoisson, link: log Response: ocorr Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 23 28.663 grupo 3 11.774 20 16.889 0.0008834 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 plot(mod2) ##Teste de comparação múltipla summary(glht(mod2, linfct=mcp(grupo="Tukey"))) Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: Tukey Contrasts Fit: glm(formula = ocorr ~ grupo, family = "quasipoisson", data = trat1_90antes) Linear Hypotheses: Estimate Std. Error z value Pr(>|z|) humi - ceifa == 0 -1.0986 0.3249 -3.382 0.00405 ** inter - ceifa == 0 -0.8979 0.3019 -2.974 0.01491 * xilo - ceifa == 0 -0.6568 0.2780 -2.363 0.08255 . inter - humi == 0 0.2007 0.3794 0.529 0.95122 xilo - humi == 0 0.4418 0.3606 1.225 0.60632 xilo - inter == 0 0.2412 0.3401 0.709 0.89176 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Adjusted p values reported -- single-step method) > ##Resultados > tapply(trat1_90antes$ocorr, trat1_90antes$grupo,mean, na.rm = TRUE)##Media dos tratamentos ceifa humi inter xilo 4.500000 1.500000 1.833333 2.333333 > tapply(trat1_90antes$ocorr,trat1_90antes$grupo,sd,na.rm=TRUE)/sqrt(tapply(trat1_90antes$ocorr,trat1_90antes$grupo,length))###Erros padrão da media ceifa humi inter xilo 0.6708204 0.5000000 0.5426274 0.3333333 *Aqui tornou-se difícil saber a diferença entre os grupos funcionais em termos de atribui a, b, c, as médias, então resolvi fazer uma análise gráfica usando o pacote HH:* > ##Visualizando as comparações > contrast.1<- + if.R(r=glht(mod2, linfct=mcp(grupo="Tukey")), + s=multicomp(mod2, plot=FALSE)) > plot(contrast.1) *Que me deu um gráfico pouco informativo, porém seguindo o help do mmc {HH}, encontrei que o gráfico de comparações deveria ser feito usando as funções glht.mmc e multicomp.mmc, que requerem o uso de aov(), então fiz:* > > ## Contraste por mmc > contrast.aov <- aov(ocorr ~ grupo, data=trat1_90antes) > contrast.mmc <- + if.R(r=glht.mmc(contrast.aov, linfct=mcp(grupo="Tukey")), + s=multicomp.mmc(contrast.aov, plot=FALSE)) > plot(contrast.mmc) > # *Pois bem, agora tenho a comparação entre os grupos, mas não utilizando a distribuição de erros quasipoisson e não sei o que fazer, pois o pacote HH não me permite utilizar os objetos criados usando glm(), para se fazer as comparações múltiplas. Alguém do grupo teria alguma sugestão, Obrigado, Alexandre* -- Alexandre dos Santos Engenheiro Florestal, Dr. Universidade Federal de Lavras Departamento de Entomologia Laboratório de Entomologia Florestal Caixa Postal 3037 37200-000 - Lavras/MG Fone: +55 (35) 9223-0304

2 2

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 09 Mai '12

09 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