- R-br - DINF UFPR

Re: [R-br] Como estimar os standard deviations de estimativas de herdabilidade em um teste
by Paulo Ricardo Gherardi Hein 22 Mar '11

22 Mar '11

Olá pessoal, Obrigado pelas dicas (função lmer do "lme4" e "SEI.GE"), mas o meu problema ainda não foi solucionado. Calcular o sd de herbabilidades estimadas a partir de modelos com efeito aleatório parece não ser uma tarefa simples no R. Os editores das revistas exigem o valor H² acompanhado pelo standard deviation/intervalo de confiança. Como é que os colegas que trabalham com melhoramento apresentam as suas estimativas de herdabilidade? Meus orientadores me ajudaram a calcular isso no SAS e no ASREML para outro design experimental (factorial mating design), mas eu não tenho licença para esses softwares e pretendo migrar 100% para o R até o ano que vem. Me desculpem por insistir no assunto, mas alguém da lista já calculou essas estatísticas no R? Abraços, Paulo Ricardo Gherardi Hein PhD candidate at University of Montpellier 2 CIRAD - PERSYST Department Research unit: Production and Processing of Tropical Woods TA B-40/16 bur.145 Adress: 73 rue Jean-François Breton 34398 Montpellier Cedex 5 - France phone: +33 4 67 61 44 51 / mobile: +33 6 47 90 57 60 email: paulo.hein(a)cirad.fr / skype: paulo_hein ------------------------------ *De:* paulo_hein <paulo_hein(a)yahoo.com.br> *Para:* R_STAT(a)yahoogrupos.com.br *Enviadas:* Terça-feira, 15 de Março de 2011 13:52:47 *Assunto:* [R_STAT] Re: Como estimar os standard deviations de estimativas de herdabilidade em um teste Prezado Fábio, obrigado pela sugestão. Prezado Fernando, você tem razão. Esse script retorna o valor do desvio padrão de cada componente da variância. No entanto, ao elevarmos os valores ao quadrado, a estimativa da herdabilidade é calculada apartir das variâncias, não é? O scrip retorna os valores de StdDev do efeito clone, interação clone x site e erro. Aqui vai um exmplo: 2581 > Resume1 2582 [[1]] 2583 Linear mixed-effects model fit by REML 2584 Data: Data 2585 Log-restricted-likelihood: -521.7632 2586 Fixed: Y ~ 1 + Site 2587 (Intercept) Site302 Site303 2588 60.79 -3.38 6.52 2590 Random effects: 2591 Formula: ~1 | Clone 2592 (Intercept) 2593 StdDev: *2.348988* 2595 Formula: ~1 | Site %in% Clone 2596 (Intercept) Residual 2597 StdDev: *5.427933* *7.052537* 2599 Number of Observations: 150 2600 Number of Groups: 2601 Clone Site %in% Clone 2602 10 30 Assim, estou calculando a broad-sense heritability da propriedade [[1]] como: H² = 2.348988^2 / (2.348988^2 + 5.427933^2 + 7.052537^2) O que eu preciso é calcular o standard deviation de cada estimativa de herdabilidade para apresentar os resultados como exigem as principais revistas da área. Algo do tipo: H² [[1]] = 0.065 (+ ou - SE) Meu orientador já calculou os SE usando o software ASREML, mas eu quero fazer isso no R. Abraços, Paulo Ricardo Gherardi Hein PhD candidate at University of Montpellier 2 CIRAD - PERSYST Department Research unit: Production and Processing of Tropical Woods TA B-40/16 bur.145 Adress: 73 rue Jean-François Breton 34398 Montpellier Cedex 5 - France phone: +33 4 67 61 44 51 / mobile: +33 6 47 90 57 60 email: paulo.hein@... / skype: paulo_hein > > > > > > ------------------------------ > > *De:* paulo_hein paulo_hein@... > > *Para:* R_STAT(a)yahoogrupos.com.br > > *Enviadas:* Segunda-feira, 14 de Março de 2011 5:24:20 > > *Assunto:* [R_STAT] Como estimar os standard deviations de estimativas de > > herdabilidade em um teste > > > > > > > > Olá, > > Trabalho com melhoramento de árvores (Eucalyptus) para producão de madeira. > > Atualmente estamos estimando os componentes de variância e herdabilidade de > > algumas propriedades da madeira proveniente de um teste clonal. O design > > experimental é composto por 10 clones, plantados em 3 sites. Cada clone é > > representado por 5 individuos, totalizando uma amostragem de 150 árvores (10 > > clones x 5 repetições x 3 sites). > > > > A maioria das estimativas do standard deviation de herdabilidades são > > feitas com base nas séries de Taylor, mas não sei como calculá-las. Alguém > > trabalha com esse assunto e conhece alguma solução? > > > > Para estimar os componentes de variância das 120 variáveis (estimados por > > espectroscopia de NIR) estamos usando dois modelos (com e sem interação). > > Para selecionar os modelos, estamos nos baseado no valor de BIC > > (selecionamos os menor BIC). Os cálculos são feitos a partir dos seguintes > > scripts: > > > > Data <- read.table("CNB150predNIR.csv", dec=".", sep=",", header=T) > > Data$Site <- as.factor(Data$Site) > > Data$Clone <- as.factor(Data$Clone) > > > > ## Estimating variance components > > Resume <- NULL > > Resume1 <- NULL > > for (i in 4:124){ > > Y <- Data[,i] > > lme <- lme(Y~ 1 + Site, random=~1|Clone, data=Data, method="REML") > > lme1 <- lme(Y~ 1 + Site, random=~1|Clone/Site, data=Data, method="REML") > > Resume <- c(Resume,list(lme)) > > Resume1 <- c(Resume1,list(lme1))} > > Resume > > Resume1 > > > > ## For calculating BIC > > Selection <- NULL > > Selection1 <- NULL > > for (i in 4:124){ > > Y <- Data[,i] > > lme <- lme(Y~ 1 + Site, random=~1|Clone, data=Data, method="ML") > > lme1 <- lme(Y~ 1 + Site, random=~1|Clone/Site, data=Data,method="ML") > > Selection <- c(Selection,list(lme)) > > Selection1 <- c(Selection1,list(lme1))} > > > > ## To display BIC results > > Resultat <- NULL > > for (i in 1:121){ lme0 <- Selection[[i]] > > lme1 <- Selection1[[i]] > > Resultat <- c(Resultat,list(anova(lme0,lme1)))} > > Resultat > > > > e as herdabilidades são estimadas da seguinte forma: > > > > H² (mle) = sdev(Clone)^2 / [sdev(Clone)^2 + sdev(Error)^2] ou > > H² (lme1) = sdev(Clone)^2 / [sdev(Clone)^2 + sdev(interaction)^2 + > > sdev(Error)^2] > > > > A maioria das estimativas do standard deviations de herdabilidades são > > feitas com base nas séries de Taylor, mas não sei como calculá-las. Alguém > > trabalha com esse assunto e conhece alguma solução? > > > > Atenciosamente, > > > > Paulo Ricardo Gherardi Hein > > PhD candidate at University of Montpellier 2 > > CIRAD - PERSYST Department > > Research unit: Production and Processing of Tropical Woods TA B-40/16 > > bur.145 > > Adress: 73 rue Jean-François Breton 34398 Montpellier Cedex 5 - France > > phone: +33 4 67 61 44 51 / mobile: +33 6 47 90 57 60 > > email: paulo.hein@... / skype: paulo_hein > > > > > > > > > > > __._,_.___ <ivanalaman@yahoo.com.br?subject=Res%3A%20%5BR_STAT%5D%20Re%3A%20Como%20estimar%20os%20standard%20deviations%20de%20estimativas%20de%20herdabilidade%20em%20um%20teste>| através de email<R_STAT@yahoogrupos.com.br?subject=Res%3A%20%5BR_STAT%5D%20Re%3A%20Como%20estimar%20os%20standard%20deviations%20de%20estimativas%20de%20herdabilidade%20em%20um%20teste>| Responder através da web<http://br.groups.yahoo.com/group/R_STAT/post;_ylc=X3oDMTJycDAydGVhBF9TAzk3N…>| Adicionar um novo tópico<http://br.groups.yahoo.com/group/R_STAT/post;_ylc=X3oDMTJmZGxhb2hzBF9TAzk3N…> Mensagens neste tópico<http://br.groups.yahoo.com/group/R_STAT/message/15816;_ylc=X3oDMTM3ajNyaGtk…>( 5) Atividade nos últimos dias: - Novos usuários<http://br.groups.yahoo.com/group/R_STAT/members;_ylc=X3oDMTJnZ3A1MjJ1BF9TAz…> 9 Visite seu Grupo<http://br.groups.yahoo.com/group/R_STAT;_ylc=X3oDMTJmZGEzb2RnBF9TAzk3NDkwND…> Você pode ser quem está faltando para alguém.<http://global.ard.yahoo.com/SIG=15ma38fhn/M=758712.14532720.14422460.129601…> [image: Yahoo! Grupos]<http://br.groups.yahoo.com/;_ylc=X3oDMTJlbTY0dXIzBF9TAzk3NDkwNDM1BGdycElkAz…> Trocar para: SÃ³ Texto<R_STAT-traditional@yahoogrupos.com.br?subject=Mudar+Formato+de+Envio:+Tradicional>, Resenha DiÃ¡ria<R_STAT-digest@yahoogrupos.com.br?subject=Envio+de+email:+Resenha>• Sair do grupo <R_STAT-unsubscribe(a)yahoogrupos.com.br?subject=Sair+do+grupo> • Termos de uso <http://br.yahoo.com/info/utos.html> . __,_._,___

3 2

Criar um ambiente LaTeX que possui lista no sumário
by Walmes Zeviani 22 Mar '11

22 Mar '11

Saudações, Fui questionado por um aluno se é possível criar um novo ambiente no LaTeX para inserir fotos e ter a lista de fotografias assim como temos a lista de figuras e lista de tabelas. Alguém pode citar algum pacote? Desde já grato e desculpem pela dúvida fora do escopo da [R-br]. 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 ==========================================================================

4 3

Enc: [r-br] Função para quebra de linha em título
by Fabio Mathias Corrêa 22 Mar '11

22 Mar '11

----- Mensagem encaminhada ---- De: Henrique Dallazuanna <wwwhsd(a)gmail.com> Para: R_STAT(a)yahoogrupos.com.br Enviadas: Sábado, 25 de Outubro de 2008 20:54:23 Assunto: Re: [R_STAT] Função para quebra de linha em título Marcio, uma alternativa talvez seja: nm <- c("Aaaaaa aaaaaaa aa", "Bbbbbb bbbbbb", "cccc") y <- c(10,4,9) barplot(y, names.arg = nm) strP <- .5 nm[strwidth( nm) > strP] <- gsub("( \\w*$)", "\n\\1", nm[strwidth( nm) > strP]) 2008/10/25 Marcio Nicolau <marcio.nicolau@ gmail.com> Prezados, > >A algum tempo, questionei na lista sobre a existência de uma função >para quebra de linha >automática para títulos de gráficos. > >Após a sugestão do Benilton (função gsub), segue abaixo o código... >caso necessitem. > > >## Pre-Processa titulos > >pretitulo <- function(orig, cut=55) { > nparts = ceiling(nchar( orig) / cut) > cuts = c(1,1:nparts * cut) > spt = strsplit(orig,"\\s")[[1]] > cs = cumsum(nchar( spt)) > > ids <- matrix(FALSE, nparts,length( nchar(spt) )) > for(i in 2:length(cuts) ) { > ids[i-1,] <- (cs >= cuts[i-1] & cs < cuts[i]) > } > > ch <- vector("character", nparts) > for(i in 1:nparts) { > ch[[i]] <- paste(spt[ids[ i,]],collapse = " ") > } > pretitulo <- paste(ch,collapse="\n") >} > >Não é o melhor, computacionalmente falando, mas funciona :) > >Abraços > >Marcio Nicolau >marcio.nicolau@ gmail.com > > > > >------------ --------- --------- ------ > >Links do Yahoo! Grupos > ><*> Para visitar o site do seu grupo na web, acesse: > http://br.groups. yahoo.com/ group/R_STAT/ > ><*> Para sair deste grupo, envie um e-mail para: > R_STAT-unsubscribe@ yahoogrupos. com.br > ><*> O uso que você faz do Yahoo! Grupos está sujeito aos: > http://br.yahoo. com/info/ utos.html > > > -- Henrique Dallazuanna Curitiba-Paraná -Brasil 25° 25' 40" S 49° 16' 22" O __._,_.___ Mensagens neste tópico (2) Responder (através da web) | Adicionar um novo tópico Mensagens | Arquivos | Fotos | Links | Banco de dados | Enquetes Alterar configurações via web (Requer Yahoo! ID) Alterar configurações via e-mail: Alterar recebimento para lista diária de mensagens | Alterar formato para o tradicional Visite seu Grupo | Termos de uso do Yahoo! Grupos | Sair do grupo Atividade nos últimos dias * 8 Novos usuáriosVisite seu Grupo Yahoo! Mail Conecte-se ao mundo Proteção anti-spam Muito mais espaço Yahoo! Barra Instale grátis Buscar sites na web Checar seus e-mails . Yahoo! Grupos Crie seu próprio grupo A melhor forma de comunicação . __,_._,___

4 3

Função para quebra de linha em título
by Fabio Mathias Corrêa 22 Mar '11

22 Mar '11

Fábio Mathias Corrêa Departamento de Estatística Universidade Estadual de Santa Cruz Tel.: 73-3680-5076 Cel.: 73-9991-8155 ----- Mensagem encaminhada ---- De: Marcio Nicolau <marcio.nicolau(a)gmail.com> Para: R_STAT(a)yahoogrupos.com.br Enviadas: Sábado, 25 de Outubro de 2008 19:47:08 Assunto: [R_STAT] Função para quebra de linha em título Prezados, A algum tempo, questionei na lista sobre a existência de uma função para quebra de linha automática para títulos de gráficos. Após a sugestão do Benilton (função gsub), segue abaixo o código... caso necessitem. ## Pre-Processa titulos pretitulo <- function(orig, cut=55) { nparts = ceiling(nchar( orig) / cut) cuts = c(1,1:nparts * cut) spt = strsplit(orig, "\\s")[[1] ] cs = cumsum(nchar( spt)) ids <- matrix(FALSE, nparts,length( nchar(spt) )) for(i in 2:length(cuts) ) { ids[i-1,] <- (cs >= cuts[i-1] & cs < cuts[i]) } ch <- vector("character" , nparts) for(i in 1:nparts) { ch[[i]] <- paste(spt[ids[ i,]],collapse = " ") } pretitulo <- paste(ch,collapse= "\n") } Não é o melhor, computacionalmente falando, mas funciona :) Abraços Marcio Nicolau marcio.nicolau@ gmail.com __._,_.___ Mensagens neste tópico (1) Responder (através da web) | Adicionar um novo tópico Mensagens | Arquivos | Fotos | Links | Banco de dados | Enquetes Alterar configurações via web (Requer Yahoo! ID) Alterar configurações via e-mail: Alterar recebimento para lista diária de mensagens | Alterar formato para o tradicional Visite seu Grupo | Termos de uso do Yahoo! Grupos | Sair do grupo Atividade nos últimos dias * 8 Novos usuáriosVisite seu Grupo Yahoo! Mail Conecte-se ao mundo Proteção anti-spam Muito mais espaço Yahoo! Barra Instale grátis Buscar sites na web Checar seus e-mails . Yahoo! Grupos Crie seu próprio grupo A melhor forma de comunicação . __,_._,___

6 6

Transformando dados com 2 fatores em uma tabela
by Augusto Ribas 22 Mar '11

22 Mar '11

Oi galera. Minha primeira postagem na nova lista que emoção :) Gostaria de ajuda com um problema simples que estou tendo aqui. #Eu queria fazer um plot assim: barplot(VADeaths, beside = TRUE, col = c("lightblue", "mistyrose", "lightcyan", "lavender", "cornsilk"), legend = rownames(VADeaths), ylim = c(0, 100)) title(main = "Death Rates in Virginia", font.main = 4) #Mas meu dados estão na forma de um data frame dessa forma... exemplo<-dados<-data.frame(cbind( Tratamento=rep(c("Controle","Ac10","Ac20","Estresse"),each=3) ,Placa=rep(1:3,4) ,Germinacao=round(rnorm(12,5,3),2) )) #Tem alguma maneira simples de fazer meu dataframe de exemplo virar uma tabela como a do "VADeaths" que é usado no help do barplot #pra fazer um plot daquele jeito, separando grupos de barrinhas tanto por tratamento como placa, no #meu caso, ou então fazer um barplot do jeito que quero(o do "VADeaths") direto com um dataframe assim do meu jeito? Obrigado -- Grato Augusto C. A. Ribas

2 2

Res: [R_STAT] Regressão de proporções com erro normal ou binomial
by Walmes Zeviani 21 Mar '11

21 Mar '11

Guilherme, Na regressão logística você observa significância porque o parâmetro de dispersão não é estimado mas sim conhecido/fixado como 1 (note que é aplicado o teste z). No modelo normal, além da estimação dos efeitos/coeficientes/parâmetros, é estimado a variância do erro aleatório e os erros padrões das estimativas estão em função disso (teste t). Se você usar uma family=quasibinomial você notará diferença. Esperamos você na [R-br]. À 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

Estudo de função!
by Elton Gean Araújo Araújo 21 Mar '11

21 Mar '11

Boa tarde! Gostaria de aplicar os conceitos de crescimento, decrescimento, máximos e mínimos de funções no R, poderiam me informar o pacote para isso? Desde já agradeço pela informação e pela nova lista! Abraço. Elton.

1 0

Especificação do modelo em lm e glm
by Diego Pujoni 21 Mar '11

21 Mar '11

#Olá pessoal, #Existem formas diferentes de se especificar um modelo na função #"lm" e "glm". Mas obtive os seguintes resultados: #Criando o data frame# X1=c(rep("A",4),rep("B",4)) X2=c(rep("C",2),rep("D",2),rep("C",2),rep("D",2)) Y=c(3,5,5,7,11,13,13,15) Table=as.data.frame(cbind(X1,X2,as.numeric(Y))) Table #Os quatro modos de se especificar o mesmo modelo fornecem a mesma deviance,# #ou seja, são iguais# glm(Y ~ X1+X2+X1:X2, data=Table, family=poisson())$deviance glm(Y ~ (X1+X2)^2, data=Table, family=poisson())$deviance glm(Y ~ X1*X2, data=Table, family=poisson())$deviance glm(Y ~ X1:X2, data=Table, family=poisson())$deviance #Mas há um deslocamento de -1.25 nos coeficientes no último modelo round(glm(Y ~ X1+X2+X1:X2, data=Table,family=poisson())$coefficients,digits=2) round(glm(Y ~ (X1+X2)^2, data=Table,family=poisson())$coefficients,digits=2) round(glm(Y ~ X1*X2, data=Table,family=poisson())$coefficients,digits=2) round(glm(Y ~ X1:X2, data=Table,family=poisson())$coefficients,digits=2) #O coeficiente de X1B:X2D sai como NA# #O quarto modo de especificar o modelo está errado? #Por que isso acontece? #Nos modelos abaixo: #Incluindo mas uma variável preditora em Table# X3=c("E","F","E","F","E","F","E","F") Table.X3=as.data.frame(cbind(X1,X2,X3,as.numeric(Y))) Table.X3 #As deviances e os coeficientes X1A:X3F e X1A:X3E também são iguais #além disso alguns coeficientes são deslocados de -1.4 (X1A:X2C,X1B:X2C,X1A:X1D), #mas os outros coeficientes não são iguais# glm(Y~X1+X2+X3+X1:X2+X1:X3,data=Table.X3,family=poisson())$deviance glm(Y~X1:X2+X1:X3,data=Table.X3,family=poisson())$deviance round(glm(Y~X1+X2+X3+X1:X2+X1:X3,data=Table.X3,family=poisson())$coefficients,digits=2) round(glm(Y~X1:X2+X1:X3,data=Table.X3,family=poisson())$coefficients,digits=2) #Especificando os coeficientes do primeiro modelo # Intercept = 1.16 # X1/X2| C D | X1/X3| E F | # A | 0 0.41| A | 0 0.41| # B |1.24 -0.25| B |1.24 -0.25| #Especificando os coeficientes do segundo modelo # Intercept = 2.56 # X1/X2| C D | X1/X3| E F | # A |-1.40 -0.99| A | 0 0.41| # B |-0.15 NA | B | 0 0.15| #Os coeficientes X1B:X3E e X1B:X3F mudam. #Alguém sabe aonde se encontra o meu erro? #Será que sou obrigado a incluir o termo de interação X2:X3? #Mas e se essas duas variáveis forem independentes? #Muito Obrigado # Diego PJ

2 1

Solução de Equação não Linear
by Carolina Pena 21 Mar '11

21 Mar '11

Boa tarde! Alguém sabe como consigo resolver uma equação não linear no R ? Por exemplo suponha que eu tenho a seguinte equação: y = 2,5*x + 3,6*x^2 + 7,8*x^4+ 10,0*x^12 Para um dado valor de y (exemplo y = 5 ) qual é o valor de x ? Obrigada, Carolina

4 4

Re: [R-br] [R_STAT] [Dúvida] Carregar dados sem ser na memória RAM
by Fernando Mayer 20 Mar '11

20 Mar '11

O pacote SOAR pode resolver. De uma olhada no manual e na vignette em http://cran-r.c3sl.ufpr.br/web/packages/SOAR/index.html []s, --- Fernando Mayer URL: http://sites.google.com/site/fernandomayer e-mail: fernandomayer [@] gmail.com 2011/3/20 pedrorafaelmarinho <pedrorafaelmarinho(a)yahoo.com.br>: > Pessoal estou com um problema ou talvez seja do R. Tenho que ler uma base de dados de 600 mb em formato DBF (DBase). O R pelo que ando lendo ele ler a base de dados e joga na memória RAM. Toda vez que carrego essa base de dados o pc trava a ponto de o mouse não funcionar. Existe alguma forma de fazer com que o R não carregue toda a base de dados diretamente na memória? Por exemplo a base de dados ficar no HD e somente carregar algumas variáveis a medida que vou precisando? Softwares simples como o Tabwin eu consigo acessar a base de dados e o pc não trava...O que está acontecendo? > O comando que usei para carregar a base foi read.dbf() usando o pacote foreign. > > > > ------------------------------------ > > Links do Yahoo! Grupos > > <*> Para visitar o site do seu grupo na web, acesse: > http://br.groups.yahoo.com/group/R_STAT/ > > <*> Para sair deste grupo, envie um e-mail para: > R_STAT-unsubscribe(a)yahoogrupos.com.br > > <*> O uso que você faz do Yahoo! Grupos está sujeito aos: > http://br.yahoo.com/info/utos.html > > >

3 4