- R-br - DINF UFPR

Adicionar os desvios padrão num gráfico
by Wenceslau Teixeira 12 Dez '12

12 Dez '12

Caro Walnes, Muito obrigado pela atenção. Funcionou tudo aqui ...apenas não aparece no grafico no eixo x as datas. Se voce ou outro membro da lista puder me ajudar em mais isto Segue o script que você corrigiu x <- c("30/04/2010","13/05/2010","28/05/2010", "15/06/2010","30/06/2010","15/07/2010", "30/07/2010","15/08/2010","30/08/2010", "15/09/2010","30/09/2010","15/10/2010", "30/10/2010","15/11/2010","30/11/2010", "15/12/2010","30/12/2010","15/01/2011", "30/01/2011","15/02/2011","28/02/2011", "15/03/2011","30/03/2011","15/04/2011", "30/04/2011","15/05/2011","30/05/2011", "15/06/2011","30/06/2011","15/07/2011") y <- c(0.278, 0.068, 0.064, 0.028, 0.030, 0.024, NA, 0.025, 0.039, 0.077, 0.081, NA, 0.043, NA, 0.343, NA, 0.195, 0.090, 0.027, 0.157, 0.122, 0.069, 0.197, 0.107, 0.103, 0.154,0.164, 0.051,0.243,0.145) SDFE10 <- c(0.294,0.033,0.058,0.013,0.006,0.006,NA,0.018,0.031,0,0, NA,0.021,NA,0.435,NA,0.057,0.049,0.005,0.220,0.159,0.029, 0.194,0.114,0.109,0.116,0.187,0,0.266,0.121) # colocar num data.frame, eliminar NA, converter para data da <- data.frame(x=as.POSIXlt(x, format="%d/%m/%Y"), y=y, sd=SDFE10) da <- na.omit(da) da$x plot(y~x, data=da, ylab=expression(Concentração~de~Ferro~na~solução~do~solo~(mg~L^{-1})),xlab="Dia", pch=4, type="b", xaxt="n",ylim=c(-0.1,0.8)) with(da, arrows(x, y-sd, x, y+sd, code=3, angle=90,length=0.03)) #axis.POSIXct(side=1, at=da$x, format="%d/%m") axis.POSIXct(side=1, at=seq(min(da$x), max(da$x), by="month"), format="%d/%m") Cordialmente, Wenceslau Teixeira Embrapa Solos -- Wenceslau Teixeira wgt007(a)gmail.com skype: wgt007

1 0

Forma "correta" e/ou mais "elegante"...
by Mauricio Morgado Oliveira 12 Dez '12

12 Dez '12

Gostaria de alguma dica (sugestão) para escrever de forma "correta" e/ou mais "elegante". Segue abaixo minha antiga programação no sas e no R (tentei anteriomente dentro de um for (loop), mas fica muito demorado). Att, Mauricio M. Oliveira Programação no SAS: if TP_PRODUTO='RM' then pai=COMPRESS('RM'||Cod_Rebanho||racapai||ano||nomepai,' '); else if livropai^=' ' & registropa^=' ' & registropa^='0' & registropa^='-' then pai=COMPRESS(livropai||'.'||registropa,' '); else if nomepai=' ' & tatuagempa=' ' then do; pai=COMPRESS(Cod_Rebanho||raca||ano||"DESC",' '); TP_PRODUTO='RM'; end; else pai=COMPRESS(Cod_Rebanho||racapai||tatuagempa||nomepai,' '); Script no R: Base[,"PAI"][Base[,"TP_PRODUTO"] == "RM"] <- paste("RM",Base[,"Cod_Rebanho"][Base[,"TP_PRODUTO"] == "RM"], Base[,"RACAPAI"] [Base[,"TP_PRODUTO"] == "RM"], Base[,"anoN"][Base[,"TP_PRODUTO"] == "RM"], Base[,"NOMEPAI"][Base[,"TP_PRODUTO"] == "RM"], sep="") Base[,"PAI"][Base[,"TP_PRODUTO"] != "RM" & Base[,"LIVROPAI"]!="" & Base[,"REGISTROPA"]!="" & Base[,"REGISTROPA"]!="0" & Base[,"REGISTRO"]!= "-"] <- paste(Base[,"LIVROPAI"][Base[,"TP_PRODUTO"] != "RM" & Base[,"LIVROPAI"]!="" & Base[,"REGISTROPA"]!="" & Base[,"REGISTROPA"]!="0" & Base[,"REGISTRO"]!="-"], ".", Base[,"REGISTROPA"][Base[,"TP_PRODUTO"] != "RM" & Base[,"LIVROPAI"]!="" & Base[,"REGISTROPA"]!="" & Base[,"REGISTROPA"]!="0" & Base[,"REGISTRO"]!="-"], sep="") Base[,"PAI"][Base[,"TP_PRODUTO"] != "RM" & Base[,"NOMEPAI"]=="" & Base[,"TATUAGEMPA"]==""] <- paste(Base[,"Cod_Rebanho"] [Base[,"TP_PRODUTO"] != "RM" & Base[,"NOMEPAI"]=="" & Base[,"TATUAGEMPA"]==""], Base[,"RACA"][Base[,"TP_PRODUTO"] != "RM" & Base[,"NOMEPAI"]=="" & Base[,"TATUAGEMPA"]==""],"DESC", sep="")

2 2

Adicionar os desvios padrão num gráfico
by Wenceslau 12 Dez '12

12 Dez '12

Caros membros da lista, Estou fazendo um gráfico simples (datas vs teores de ferro). Eu queria adicionar os desvios padrão das avaliações dos teores em cada data. Que foram calculados separadamente. Mas não estou conseguindo. Talvez alguém tenha experiência com estes gráficos e possa me ajudar. Segue o script ########################### # SCRIPT PARA GRAFICOS DE ELEMENTOS NO SOLO #DATAS DE COLETAS x<-1:30 x<-factor x<- c("30/04/2010","13/05/2010","28/05/2010","15/06/2010","30/06/2010", "15/07/2010", "30/07/2010" ,"15/08/2010", "30/08/2010","15/09/2010","30/09/2010","15/10/2010","30/10/2010","15/11/2010 ","30/11/2010","15/12/2010", "30/12/2010", "15/01/2011", "30/01/2011", "15/02/2011", "28/02/2011" ,"15/03/2011","30/03/2011","15/04/2011","30/04/2011","15/05/2011","30/05/201 1","15/06/2011","30/06/2011", "15/07/2011") #DADOS DE FE AVALIADOS NO LABORATORIO y<-1:30 y <- c(0.278, 0.068, 0.064, 0.028, 0.030, 0.024, NA, 0.025, 0.039, 0.077, 0.081, NA, 0.043, NA, 0.343, NA, 0.195, 0.090, 0.027, 0.157, 0.122, 0.069, 0.197, 0.107, 0.103, 0.154,0.164, 0.051,0.243,0.145) # DESVIOS PADRAO DOS TEORES DE FE SDFE10 <- c(0.294, 0.033, 0.058,0.013,0.006,0.006,NA,0.018,0.031,0,0,NA,0.021,NA,0.435,NA,0.057,0.049, 0.005,0.220,0.159,0.029,0.194,0.114,0.109,0.116,0.187,0,0.266,0.121) #PLOTAR O GRAFICO plot(y, ylim=c(-0.1,0.7), axes=F, ylab="Concentração de Ferro na solução do solo (mg L-1)", xlab="Dia", pch=4, type="b") #Adicionar os desvios padrão (SDFE10) arrows(y,y-SDFE10,y, y+SDFE10, code = 3, angle = 90,length = 0.03) Wenceslau Geraldes Teixeira Embrapa Solos (www.cnps.embrapa.br) Rua Jardim Botânico, 1024 CEP 22460-000 Rio de Janeiro - RJ - Brasil Tel: 55 21 21794583 Skype: wgt007 E-mail: wenceslau(a)cnps.embrapa.br / wenceslau.teixeira(a)pq.cnpq.br CV: http://lattes.cnpq.br/4136381865367906 ________________________________ II Encontro Brasileiro de Física do Solo 2sd Brazilian Soil Physics Meeting Rio de Janeiro - 06 a 10 de Maio de 2013 http://www.bspm2013.org/

2 1

Imputação múltipla: como fazer operações com os dados imputados.
by Lucas Petri Damiani 12 Dez '12

12 Dez '12

Boa noite, sou leitor assíduo da lista, embora nunca tenha participado das discussões. Porém, hoje vou me arriscar e passar um problema pra vcs. Espero que possam me ajudar. Estou realizando uma análise num banco de dados de cerca de 10 mil pacientes de Trauma. O objetivo da pesquisadora é testar a validade de um escore de gravidade (preditor de sobrevida) em pacientes que entram no pronto socorro e eventualmente promover alterações nesse escore. Porém, algumas medidas que ela deseja testar possuem cerca de 20% de missing. E a exclusão desses pacientes faz com que o novo escore predito fique um bocado diferente do original, o que me fez ir atrás de alternativas para tratar esses dados faltantes e assim me arrisquei nesse pacote *mice. *Em testes mais simples que fiz, fiquei bastante empolgado com o método. Enfim, tentei reproduzir abaixo o problema que encontrei. ##### Banco de dados de exemplo: dados <- matrix(c( 4, 6, 5, 34, 120, 5, 1, 0, NA, 1, 1, 3, 80, 6, 0, 1, 2, 4, 2, 24, 90, 5, 1, 1, 3, 5, 5, 5, 98, 5, 1, 1, NA, 6, 5, 40, 94, 4, 0, 0 , NA, 6, 5, 12, 87, 3, 0, 1, 1, 5, NA, 32, 70, 2, 0, 1, 1, 5, 1, 7, 60, 3, 1, 0, 2, NA, 2, 25, 80, 5, 1, 0, 4, 6, 5, 24, 101, 7, 1, 0, 1, 5, 1, 1, 120, 8, 1, 0, 3, 5, 3, 12, 200, 6, 0, 1, 3, 3, 1, 3, 40, 5, 1, 0, NA, 5, 1, 35, 60, 4, 0, 1, 1, 4, 1, 32, 87, 3, 1, 0, 4, 5, 1, 23, 30, 2, 1, 1, 2, 1, 1, 5, 46, 1, 0,1), byrow=T, ncol=8) dados <- as.data.frame(dados) colnames(dados) <- c("ECGL_AO", "ECGL_MRM", "ECGL_MRV", "FR", "PAS", "GRAU", "IOT", "SOBREVIDA") dados$IOT <- as.factor(dados$IOT) dados$ECGL <- apply(dados[, c("ECGL_AO", "ECGL_MRV", "ECGL_MRM")],1,sum) dados$ECGL_cut <- cut(dados$ECGL, breaks=c(3,4,6,9,13,16), right=FALSE) # Notem que a medida RTS é um escore derivado de outras medidas, que pretendo imputar. dados$RTS <- (0.9368*(dados$ECGL_cut=="[4,6)") + 2*0.9368*(dados$ECGL_cut=="[6,9)") + 3*0.9368*(dados$ECGL_cut=="[9,13)") + 4*0.9368*(dados$ECGL_cut=="[13,16)") + 0.7326*(dados$PAS>=1 & dados$PAS<=49) + 2*0.7326*(dados$PAS>=50 & dados$PAS<=75) + 3*0.7326*(dados$PAS>=76 & dados$PAS<=89) + 4*0.7326*(dados$PAS>=90) + 0.2908*(dados$FR>=1 & dados$FR<=5) + 2*0.2908*(dados$FR>=6 & dados$FR<=9) + 4*0.2908*(dados$FR>=10 & dados$FR<=29) + 3*0.2908*(dados$FR>=30)) #Imputando dados nas variáveis "ECLG": require(mice) PRED <- (1 - diag(1, ncol(dados[,1:8]))) PRED[8,1:8] <- PRED[1:8,8] <- 0 teste <- mice(dados[,1:8], method=c("pmm","pmm", "pmm", "", "", "", "", ""), seed=15, predictorMatrix = PRED) #O meu problema começa agora. Eu não sei como calcular medidas derivadas no objeto teste. E a fórmula abaixo retorna erro, já conferi várias vezes muito dos parênteses e etc, mas creio que o problema seja que a fórmula é muito grande entro da função I(x), pois quando eu reduzo o modelo, eu consigo resultados satisfatórios. Me corrijam se estiver errado, por favor. fit_error <- with(teste, glm(SOBREVIDA ~ I(0.9368*((ECGL_AO + ECGL_MRV + ECGL_MRM) >= 4) & ((ECGL_AO + ECGL_MRV + ECGL_MRM) <= 5) + 2*0.9368*(((ECGL_AO + ECGL_MRV + ECGL_MRM) >= 6) & ((ECGL_AO + ECGL_MRV + ECGL_MRM) <= 8)) + 3*0.9368*I(((ECGL_AO + ECGL_MRV + ECGL_MRM) >= 9) & ((ECGL_AO + ECGL_MRV + ECGL_MRM) <= 12)) + 4*0.9368*I(((ECGL_AO + ECGL_MRV + ECGL_MRM) >= 13) & ((ECGL_AO + ECGL_MRV + ECGL_MRM) <= 15)) + 0.7326*I(PAS>=1 & PAS<=49) + 2*0.7326*I(PAS>=50 & PAS<=75) + 3*0.7326*I(PAS>=76 & PAS<=89) + 4*0.7326*I(PAS>=90) + 0.2908*I(FR>=1 & FR<=5) + 2*0.2908*I(FR>=6 & FR<=9) + 4*0.2908*I(FR>=10 & FR<=29) + 3*0.2908*I(FR>=30)) + IOT, family="binomial")) #Enfim, gostaria de rodar algo como isso escrito abaixo, mas não sei como definir essa variável RTS dentro do arquivo teste. fit <- with(teste, glm(SOBREVIDA ~ RTS + IOT, family="binomial")) fit2 <- glm(SOBREVIDA ~ RTS + IOT, family="binomial", data=dados) pool(fit) # Uma saída que imagino é gerar os cinco bancos imputados através da função "complete(teste, m)" e ajustar a regressão. Mas dessa forma não posso utilizar a função *pool. *Teria que reprogramá-la. De qualquer forma, estou perdido. Qualquer ideia é bem vinda. Abraços, Lucas Petri Damiani

2 2

Re: [R-br] Trabalhar com model.matrix
by Walmes Zeviani 11 Dez '12

11 Dez '12

A model.matrix() te dá a matriz sob a restrição paramétrica considerada no ajuste, que por padrão é a contraste tratamento (anula 1 nível). Para ter a matriz sem restrições você vai ter que fazer a matriz para cada termo do seu modelo, removendo o intercepto e colando lado à lado, por exemplo mA <- model.matrix(~-1+A) mB <- model.matrix(~-1+B) mAB <- model.matrix(~-1+A:B) m <- cbind(A, B, AB) Submeta suas dúvidas para lista de discussão de R, mais pessoas podem contribuir, podem haver diferentes soluções, respostas em menos tempo, as mensagens ficam arquivadas para referência/consulta futura, as mensagens podem ser úteis para outros usuários. 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 skype: walmeszeviani twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================

1 0

Teste no R
by Márcia Varela da Silva 11 Dez '12

11 Dez '12

Não consigo + criar um novo tópico na R-br. Qual comando para aplicação do Teste de Bonferroni no R (para MANOVA - Comparações múltiplas multivariadas)? Desde já agradeço Att. Márcia Varela da SilvaMestranda em Produção Vegetal (Fitossanidade) - Universidade Federal do Espírito Santo(28)9964-3488(27)9812-1832

2 1

Grupo de estudo de Geoestatistica
by Hélio Gallo Rocha 11 Dez '12

11 Dez '12

Prezados, Não sei é prática desta lista, mas haveriam interessados em formar um grupo de estudo especifico em geoestatística? Hélio

15 17

Alterar source code
by Davi 10 Dez '12

10 Dez '12

Boa tarde! Estou precisando alterar o source code da função readline(), do base package. A questão é que preciso de mais de 256 caracteres e nos detalhes da ajuda da função (abaixo) parece que isto é possível de ser feito. Alguém sabe como? Obrigado! Details The prompt string will be truncated to a maximum allowed length, normally 256 chars (but can be changed in the source code). ___________________________________________________ Davi Butturi-Gomes Doutorando em Estatística e Experimentação Agronômica - ESALQ /USP - Piracicaba Mestre em Biometria - IB/UNESP - Botucatu Ecólogo - IB/UNESP - Rio Claro Skype: davibg87 (19)8206.6565

2 1

Re: [R-br] Shapefile
by Thiago V. dos Santos 10 Dez '12

10 Dez '12

Taynãna, Estou respondendo com cópia para a lista, pois esse assunto tem despertado interesse de muitos usuários que podem se beneficiar com a solução da sua dúvida. Não é possível entender completamente o seu problema ser ter amostra dos dados. Se você tem um objeto da classe SpatialPolygon que precisa ser convertido para SpatialPolygonDataFrame para ser escrito em shapefile, você pode usar a função 'SpatialPolygonsDataFrame' do pacote sp. Você precisa conhecer a natureza dos seus dados para saber o que deve fazer. Vou postar abaixo um exemplo de conversão entre classes SpatialPolygon* que encontrei no manual do pacote maptools. Esse exemplo deve te ajudar: library(maptools) xx <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1], IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66")) SPxx <- as(xx, "SpatialPolygons") class(SPxx) IDs <- sapply(slot(SPxx, "polygons"), function(x) slot(x, "ID")) df <- data.frame(rep(0, length(IDs)), row.names=IDs) SPDFxx <- SpatialPolygonsDataFrame(SPxx, df) class(SPDFxx) summary(SPDFxx) tf <- tempfile() writePolyShape(SPDFxx, tf) getinfo.shape(tf) Espero que ajude, Thiago. ________________________________ From: Taynãna César Simões <tay_estatistica(a)yahoo.com.br> To: Thiago V. dos Santos <thi_veloso(a)yahoo.com.br> Sent: Monday, December 10, 2012 8:59 AM Subject: Re: [R-br] Shapefile Thiago e Paulo, Ainda não consegui transformar o objeto SpatialPolygons em um shapefile. Podem me dar um exemplo de como colocar os argumentos? O help não ajudou muito. Obrigada, Taynãna ________________________________ De: Thiago V. dos Santos <thi_veloso(a)yahoo.com.br> Para: "r-br(a)listas.c3sl.ufpr.br" <r-br(a)listas.c3sl.ufpr.br>; Taynãna César Simões <tay_estatistica(a)yahoo.com.br> Enviadas: Sexta-feira, 7 de Dezembro de 2012 15:21 Assunto: Re: [R-br] Shapefile Taynãna, Há funções apropriadas para escrever shapefiles nos pacotes maptools(writePolyShape) e rgdal (writeOGR), mas você provavelmente terá que converter seus dados em SpatialPolygonsDataFrame antes. Thiago. ________________________________ From: Taynãna César Simões <tay_estatistica(a)yahoo.com.br> To: "r-br(a)listas.c3sl.ufpr.br" <r-br(a)listas.c3sl.ufpr.br> Sent: Friday, December 7, 2012 2:45 PM Subject: [R-br] Shapefile Pessoal, Alguém sabe como exportar um objeto da classe SpatialPolygons como um arquivo shapefile? Att, Taynãna Fiocruz/RJ _______________________________________________ R-br mailing list R-br(a)listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível. _______________________________________________ R-br mailing list R-br(a)listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.

1 0

Implementação Função
by FHRB Toledo 10 Dez '12

10 Dez '12

Senhores, Dia ou outro me vejo com a necessidade de aplicar o teste de Bartlett (função bartllet.test(), do pacote stats). Em um desses usos regulares notei certo "comportamento" anômalo na mesma. O uso que dou a ela era para testar a homogeneidade de variâncias entre modelos ajustados, ou seja, aplicava-a sobre uma lista de saídas da função lm(), da forma: bartlett.test(list(x = lm(...), z = lm(...), w = lm(...), y = lm(...))), De acordo com a sintaxe apresentada no help() da mesma ( http://stat.ethz.ch/R-manual/R-patched/library/stats/html/bartlett.test.html ). Para investigar o porque desse comportamento anômalo, além do help() indicado no link, busquei a citação do help, de autoria próprio Bartlett ( http://www.jstor.org/stable/96803) E a fonte da função disponível no tar.gz do R 2.15.2, que baixei do CRAN ( http://cran.r-project.org/src/base/R-2/R-2.15.2.tar.gz, acesso hoje 10 de dezembro de 2012). Transcrevo no final da mensagem um trecho da fonte, que se encontra em: ~/R-2.15.2/src/library/stats/R/bartlett.test.R. Pois bem, lendo a código fonte, o artigo do Bartlett e consultando o Livro Steel, Torrie and Dickey "Principles and Procedures of Statistics: A biometrial Approach", terceira edição, capítulo 19, página 480. Pude perceber que todo o fundamento do teste aplica-se a variâncias. Entretanto ao observar a implementação da função (linha destacada) observei que ao "capturar" sum(obj$residuals^2), a função estaria atendo-se a soma de quadrados do resíduo, e não a variância. Obviamente, uma soma de quadrados, dividido pelo seu respectivo grau de liberdade é uma quadrado médio e por sua vez uma variância, entretanto, isso não é realizado! Pelos testes que fiz, ao utilizar a função com diferentes modelos ajustados com o mesmo grau de liberdade dos resíduos, pela proporcionalidade das somas de quadrados o resultado do teste é satisfatório, entretanto, para diferentes modelos, com diferentes graus de liberdades o resultado do teste não confere com o obtido calculado a mão! Diante do exposto, gostaria de ouvir a opinião dos senhores quanto a esse fato, estou certo? Em caso positivo, qual seria um meio elegante de entrar em contato com o R CORE TEAM no sentido de apresentar esse fato? Agradeço qualquer explicação, sugestão, etc. att, FH bartlett.test.default <- function(x, g, ...) { LM <- FALSE if (is.list(x)) { if (length(x) < 2L) stop("'x' must be a list with at least 2 elements") DNAME <- deparse(substitute(x)) if (all(sapply(x, function(obj) inherits(obj, "lm")))) LM <- TRUE else x <- lapply(x, function(x) x <- x[is.finite(x)]) k <- length(x) } else { if (length(x) != length(g)) stop("'x' and 'g' must have the same length") DNAME <- paste(deparse(substitute(x)), "and", deparse(substitute(g))) OK <- complete.cases(x, g) x <- x[OK] g <- factor(g[OK]) k <- nlevels(g) if (k < 2) stop("all observations are in the same group") x <- split(x, g) } if (LM) { n <- sapply(x, function(obj) obj$df.resid) v <- sapply(x, function(obj) sum(obj$residuals^2)) ***** DESTAQUE ***** DESTAQUE ***** DESTAQUE ***** } else { n <- sapply(x, "length") - 1 if (any(n <= 0)) stop("there must be at least 2 observations in each group") v <- sapply(x, "var") } n.total <- sum(n) v.total <- sum(n * v) / n.total STATISTIC <- ((n.total * log(v.total) - sum(n * log(v))) / (1 + (sum(1 / n) - 1 / n.total) / (3 * (k - 1)))) PARAMETER <- k - 1 PVAL <- pchisq(STATISTIC, PARAMETER, lower.tail = FALSE) names(STATISTIC) <- "Bartlett's K-squared" names(PARAMETER) <- "df" RVAL <- list(statistic = STATISTIC, parameter = PARAMETER, p.value = PVAL, data.name = DNAME, method = "Bartlett test of homogeneity of variances") class(RVAL) <- "htest" return(RVAL) }

1 0