- R-br - DINF UFPR

Canal #R-BR no IRC
by Jodavid Ferreira 08 Nov '13

08 Nov '13

Para aqueles que utilizam IRC, agora existe o canal *#r-br* no servidor* irc.freenode.net:6667 <http://irc.freenode.net:6667>* o canal *#R* já existe é bastante utilizado e acessado, agora existe também o canal brasileiro. -- Jodavid Ferreira Responsável pela comissão de Comunicação do MJ da RCCPB - http://rccpbjovem.com Membro da equipe Estadual do Ministério de Comunicação da RCCPB - http://rccpb.com.br Responsável pelo Ministério de Comunicação do Grupo de Oração Boa Esperança - http://grupoboaesperanca.com Bacharelando em Estatística - UFPB Estagiário no Laboratório de Estatística Aplicada ao Processamento de Imagens e Geoprocessamento - http://de.ufpb.br/~leapig Site Pessoal - http://jodavidferreira.com Twitter - @JodavidFerreira <http://twitter.com/JodavidFerreira>

5 4

Curso R
by Rodolfo Timoteo da Silva 08 Nov '13

08 Nov '13

-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA256 Apenas umas doletas: http://www.jigsawacademy.com/analytics-courses/r-analytics Att - -- - --- Rodolfo Timóteo da Silva Linux Counter: 359362 skype: zhushazang Ribeirão Preto - SP "Qualidade é inversamente proporcional a variabilidade" - Montgomery - --- -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux) iF4EAREIAAYFAlJ8zV8ACgkQSYE5fAFtI6tR3AD+OhGJIZTEYj6+XKdnGBR7fugr czqJsO3UZx1CW1k2oPwA/0yqdkWB2cPX0H+eF/QB/V10azlRJqgshfT/5pfk0+aI =CBdS -----END PGP SIGNATURE-----

1 0

Ortogonalização Gram-Smith
by Carlos Pereira 08 Nov '13

08 Nov '13

exite outra função alem de qr() para fazer a ortogonalização Gram-Smith Carlos Pereira da Silva Bacharel em Estatística – UEPB Mestrando em Estatística e Experimentação Agropecuária - UFLA Lavras - MG Cel: 035 - 9119 - 2245 Email: ccpsilva81(a)hotmail.com __________________________________________ Tenha Um Dia Abençoado Por DEUS!!! Sl – 34, 8

1 0

Monitorar pacote no CRAN.
by Pedro Rafael 07 Nov '13

07 Nov '13

Pessoal vocês conhecem alguam solução para monitorarmos nossos pacote no CRAN? Por exemplo saber quantas pessoas baixaram o pacote, localidade, etc. Os artigos que utilizam um pacote do R são obrigados a citar esse pacote? Saudações, Pedro Rafael.

4 6

[Ajuda] Pacote xlsx
by Andre Oliveira 07 Nov '13

07 Nov '13

Pessoal o que etsou fazendo de errado ao tentar carragar meus dados? O que bseria este indice? require(xlsx) dados=read.xlsx("teste.xlsx") Erro em read.xlsx("teste.xlsx") : Please provide a sheet name OR a sheet index. obrigado

9 13

Acentuação em Markdown
by Fernando Souza 07 Nov '13

07 Nov '13

Caros amigos, tisabe me informar se é possível e como se faz para colocar acentos em texto Markdown. Estou utilizando o Rstudio presentation e sempre que ustilizo acentos um sinal de interrogação aparece no lugar. Como soluciono isso? att -- ===================================================== Fernando Souza Zootecnista, DSc. Produção Animal cel: (+55) 82 8113-8781 e-mail:nandodesouza@gmail.com ==================================================

1 0

% de cada area na krigagem
by Hélio Gallo Rocha 06 Nov '13

06 Nov '13

Caros da lista, Situação: mapa de krigagem, dividido em 5 níveis ou poderia ser em quantilis Problema: Teria alguma função pronta com retorno da % de cada nível? grato -- Hélio Gallo Rocha IFSULDEMINAS - Câmpus Muzambinho

3 5

hidden function
by Fernanda De Bastiani 06 Nov '13

06 Nov '13

Caros, minha dúvida está relacionada a "hidden functions", em como acessá-las. Por exemplo: se tenho a funcão f() e eu não tenho interesse em colocar no NAMESPACE em export(f), ela será uma "hidden function", que pode ser acessada através do comando packagename:::f. No entanto, o R tem enviado alguns warnings em relação a este comando com ::: Alguém sabe me dizer a melhor forma para acessar este tipo de função? Muito obrigada, Fernanda

2 5

Selecionar posição de um número
by Jodavid Ferreira 06 Nov '13

06 Nov '13

Gostaria de saber como faço para selecionar o valor. Por Exemplo: 235478545 365984594 Gostaria de selecionar o 3 no primeiro valor e o 6 no segundo valor, ou seja, sempre a segunda posição da esquerda para a direita, como posso fazer isso? -- Jodavid Ferreira Responsável pela comissão de Comunicação do MJ da RCCPB - http://rccpbjovem.com Membro da equipe Estadual do Ministério de Comunicação da RCCPB - http://rccpb.com.br Responsável pelo Ministério de Comunicação do Grupo de Oração Boa Esperança - http://grupoboaesperanca.com Bacharelando em Estatística - UFPB Estagiário no Laboratório de Estatística Aplicada ao Processamento de Imagens e Geoprocessamento - http://de.ufpb.br/~leapig Site Pessoal - http://jodavidferreira.com Twitter - @JodavidFerreira <http://twitter.com/JodavidFerreira>

2 2

Dúvida Knitr e Markdown
by Fernando Souza 05 Nov '13

05 Nov '13

Estou fazendo uma apresentação utilizando a ferramenta presentation do Rstudio o qual utiliza Markdown e knitr. eu fiz um gráfico e gostaria de aumentar o tamnho de sua saída no slide. Alguém sabe como se faz isso? E como centralizar o texto no slide? -- ===================================================== Fernando Souza Zootecnista, DSc. Produção Animal cel: (+55) 82 8113-8781 e-mail:nandodesouza@gmail.com ==================================================

1 0

Operação com expand.grid
by ASANTOS 05 Nov '13

05 Nov '13

Boa tarde Pessoal, Estou com dificuldades de fazer um expand.grid() alterando-se as distância, como por exemplo se faço: # Ponto inicial p0<-c(1,1) c0 <- seq(p0[1],p0[1]+190,by=2) d0 <- seq(p0[2],p0[2]+47,by=3) d1 <- expand.grid(x=c0, y=d0) ## Coordenadas do grid espaçadas 3 x 2 plot(d1[,1],d1[,2]) #---------- Crio 16 pontos no eixo y e 32 pontos no eixo x, totalizando 1536 pontos, no entanto, eu gostaria de criar os mesmo 1536 pontos, mas gostaria de fazer da seguinte maneira: a cada 8 pontos na linha e 8 na coluna, distanciar em 4 unidades na linha e na coluna do próximo agrupamento de 8 pontos na linha e 8 na coluna para todos os pontos, formando ao final 24 agrupamentos de 8 pontos na linha por 8 na coluna separados de maneira equidistantes em 4 unidades. Alguém poderia me dar uma ajuda, pois fiz inúmeras tentativas sem sucesso? 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 ======================================================================

7 8

pseudo-inversa
by joão Rodrigo Castro 05 Nov '13

05 Nov '13

Pessoal, desejo fazer uma análise discriminante de Fisher e um dos passos é calcular a a inversa da matriz de covariância conjunta entre duas matrizes grandes. Essa matriz de covariância conjunta é singular e por isso não dá para calcular sua inversa, daí recebi uma orientação de usar a inversa de Moore-Penrose. Alguém sabe como efetuar esse cálculo no R? -- João Rodrigo de Castro Bolsista Programa de Educação Tutorial - PET Graduando em Meteorologia Universidade Federal de Pelotas

2 2

Dúvida Tutorial Geoestatística
by Helder Gramacho 04 Nov '13

04 Nov '13

Boa noite, sou iniciante na Geoestatística e no R, estava praticando numa espécie de tutorial situado no seguinte link: http://jurachumbo.blogspot.com.br/p/resultados_06.html Estava executando os comandos até que apareceu um erro que não estou conseguindo entender, seguem os comandos e o erro: data(jura) jura.grid=juragrid.dat coordinates(jura.grid)=~Xloc+Yloc gridded(jura.grid) = TRUE m<-vgm(0.12645599,"Exp", 0.1593586, 0.04550649) x<-krige(log(Pb)~1,jura, jura.grid, model=m) Error in krige(log(Pb) ~ 1, jura, jura.grid, model = m) : erro na avaliação do argumento 'locations' na seleção do método para a função 'krige': 'Error: object 'jura' not found ' Alguém saberia me explicar qual é o problema? Desde já agradeço, Hélder Gramacho dos Santos Engenheiro Agrônomo Mestrando Ciênc. Geod. Tec. Geoinformação-UFPE *agrohelder(a)gmail.com <agrohelder(a)hotmail.com>*

3 2

Duvida - Ellipse
by Carlos Pereira 04 Nov '13

04 Nov '13

Preciso construir elipses de confiança para um conjunto de dados bidimensionais. Mas não quero assumir normalidade. Conheço as posteriores condicionais dos parâmetros, mas meu interesse e em uma variável defendida como a multiplicação da raiz quadrada dos valores singulares pelas coordenadas dos autovetores. Se eu fizer esta multiplicação nas cadeias consigo gerar uma amostra a posteriori, mas não conheço a distribuição conjunta destas variáveis. Desta forma tenho que construir elipses baseadas em agrupamento de pontos. E não sei como construí- las. Logo se alguém souber como posso fazer peço que me ajude. Carlos Pereira da Silva Bacharel em Estatística – UEPB Mestrando em Estatística e Experimentação Agropecuária - UFLA Lavras - MGEmail: ccpsilva81(a)hotmail.com __________________________________________ Tenha Um Dia Abençoado Por DEUS!!! Sl – 34, 8

3 2

Codificação de arquivo .nnw
by Andre Oliveira 04 Nov '13

04 Nov '13

Boa noite, ha 3 anos contrui uma apostila usando o Latex + Sweave .rnw e foi criado no Windows com TinnR. Resolvi atualisar a apostila e compilar o arquivo .rnw, porém, percebi que meu arquivo .rnw está com toda ascentuação comprometida. Onde tem ascentos, til tudo com ponto de interrogação. Recorro a vocês do grupo para tentar resolver este problema. Como mudei para Linux heguei abrir este arquivo no Linux com Rstudio. Alguma dica? obrigado André Oliveira Souza

4 4

retidada de tendencia sem UTM
by Hélio Gallo Rocha 04 Nov '13

04 Nov '13

Caros da lista, Estudando outros dados, não houve problema na validação cruzada para retirada de tendencia polinomial de primeiro grau, achei que fosse pela quantidade de pontos, mais de 300. Como minhas coordenadas estão em UTM, subtrai todas elas pelo menor valor. Rodei o CRM... sem problemas. Para quem está em fase de estudos da geoestatistica, ainda faltam pontos de dúvida, favor darem suas opiniões, listo aqui minhas dúvidas e acrescentem outras: 1. Sempre temos de retirar tendencia? 2. Poderia retirar tendencia com uma covariavel fenológica da planta? penso que sim, já que a retirada polinomial de primeiro grau utiliza as coordenadas. 3. Alguém teria algum artigo para embasar este tipo de situação? Grato a todos -- Hélio Gallo Rocha IFSULDEMINAS - Câmpus Muzambinho

2 3

Inserir desvio padrão em gráfico de barras
by Rodrigo Muniz 04 Nov '13

04 Nov '13

Prezados membros da lista gostaria de pedir o auxílio de vocês. A partir dos seguintes comando: dados=cbind(c("Mata Nativa", "Citrus 37", "Citrus 37", "Citrus 24", "Citrus 18", "Citrus 31","Pasto 37"), tapply(arenosa[1:21, "ECA"], arenosa[1:21, "Trat"], mean), tapply(arenosa[1:21, "ECB"], arenosa[1:21, "Trat"], mean), id=1:nrow(tapply(arenosa[1:21, "ECA"], arenosa[1:21, "Trat"], mean))) colnames(dados)<-c("Classes", "0-10", "10-20", "id") dados= as.data.frame(dados) barchart(value~Classes, groups=variable, data=melt(dados, measure.var=2:3), auto.key=list(columns=2), horizon=FALSE, ylab="Estoque Carbono", main="Área 3") Eu gero a seguinte tabela de dados com o respectivo gráfico de barra. Classes 0-10 10-20 id 301 Mata Nativa 10.04 7.6 1 312 Citrus 37 12.36 23.12 2 322 Citrus 37 15.32 27.58 3 332 Citrus 24 4.38 8.85 4 342 Citrus 18 10.43 21.57 5 352 Citrus 31 9.19 18.03 6 364 Pasto 37 6.69 12.53 7 Eu gostaria de inserir neste gráfico de barras o desvio padrão de cada barra. Como eu poderia fazer isso? Segue a tabela com os respectivos desvios padrão. 0-10 10-20 Mata Nativa 3.4250012 0.7715569 Citrus 37 2.0617064 4.1486182 Citrus 37 2.8501287 4.7112737 Citrus 24 0.1931321 0.3666515 Citrus 18 1.7760349 4.9343726 Citrus 31 0.6730775 1.6246333 Pasto 37 2.4274060 2.8266116 Desde já agradeço. -- *Rodrigo A. Muniz* Eng. Agrônomo. Ms Produção Vegetal (UENF) Doutorando em Engenharia de Sistemas Agrícolas (ESALQ/USP) Cel (19) 8300-4333 (Pessoal) Tel (19) 3375-1196 (Residêncial) Tel (19) 3447-8537 (Trabalho)

2 1

Dúvida no Mongo DB
by Savano Pereira 04 Nov '13

04 Nov '13

Prezados, Estou acessando um MongoDB remotamente utilizando o R. Gero uma consulta, porém, quando crio um data frame com a consulta ele está vazio, alguém sabe me dizer onde está o problema?. Estou começando a me inteirar sobre o MongoDB. Envio o log do programa... Desde já agradeço. > library(rmongodb) > > # parâmetros > host <- "xxxx" > username <- "xxx" > password <- "xxxx" > db <- "censo2010" > > #conectando a base > teste <- mongo.create(db="censo2010", host="xxxx", username="xxxx",password="xxxx") > > #Obtendo a lista de bancos no MongoDBserver > print(mongo.get.databases(teste)) NULL > > #obtendo a lista de collection > print(mongo.get.database.collections(teste, db)) [1] "censo2010.censo2010" > > #Returns a fresh mongo.bson.buffer object ready to have data appended onto it > buf <- mongo.bson.buffer.create() > > #Start a subobject within a mongo.bson.buffer > mongo.bson.buffer.start.object(buf, 'V0001') [1] TRUE > #Append a name/value pair into a mongo.bson.buffer > mongo.bson.buffer.append(buf, '$eq', '11') [1] TRUE > #Finish a subobject or array within a mongo.bson.buffer > mongo.bson.buffer.finish.object(buf) [1] TRUE > > #Convert a mongo.bson.buffer object to a mongo.bson object > query <- mongo.bson.from.buffer(buf) > > #Find records in a collection > cur <- mongo.find(teste, 'censo2010.censo2010', query = query) > > V0001 <- NULL > while (mongo.cursor.next(cur)) { + value <- mongo.cursor.value(cur) + V0001 <- rbind(V0001, mongo.bson.value(value, 'V0001')) + } > data2 <- data.frame(V0001=V0001) # AQUI O DATAFRAME ESTÁ VAZIO -- Savano Pereira TEL: 21 8771-9063

1 0

R Markdown através do Rcmdr
by pjobenil＠gmail.com 04 Nov '13

04 Nov '13

Boa tarde pessoal Passei um trabalho simples de análise de dados para uns alunos de engenharia e espero que eles usem o R Markdown através do Rcmdr para gerar um relatório com as análise para me entregar. Porém o ambiente não aceita caracteres latinos (ã, õ, í, etc) e o texto em html sai com uma série de cacos. Alguém conhece alguma forma de solucionar isso? Agradeço antecipadamente, Atenciosamente Jobenil Júnior Enviado do Email do Windows

4 6

Interface em R
by Tiago Souza Marçal 04 Nov '13

04 Nov '13

Boa tarde pessoal, Andei dando uma olhada no pacote tclk e suas funcionalidades na criação de interfaces. No entanto, gostaria de saber se existe algum plugin que possa ser instalado no computador para que a interface pudesse ser construida via netbeans. Att. Tiago. ################################################################# Tiago de Souza Marçal - Graduando em Agronomia pelo CCA-UFES Bolsista de Iniciação Científica da área de Genética e Melhoramento de Plantas #################################################################

3 3

[Ajuda] Arquivo Sweave TexLive 2013
by Andre Oliveira 04 Nov '13

04 Nov '13

Pessoas, querendo atualizar no Ubtntu12.04 instalei o TexLive 2013 era TexLive 2009, porém, meu arquivo .rnw não acha a biblioteca Sweave. Achei estarnaho, pois, nunca tive problema com Sweave no Linux. Alguém sabe onde devo coloca o arquivo Sweave em que pasta do TeXLive 2013? André Oliveira Souza

1 0

Erro na função Krigagem Ordinária
by Helder Gramacho 03 Nov '13

03 Nov '13

Boa noite, sou iniciante na Geoestatística e no R, estava praticando numa espécie de tutorial situado no seguinte link: http://jurachumbo.blogspot.com.br/p/resultados_06.html Estava executando os comandos até que apareceu um erro que não estou conseguindo entender, seguem os comandos e o erro: data(jura) jura.grid=juragrid.dat coordinates(jura.grid)=~Xloc+Yloc gridded(jura.grid) = TRUE m<-vgm(0.12645599,"Exp", 0.1593586, 0.04550649) x<-krige(log(Pb)~1,jura, jura.grid, model=m) Error in krige(log(Pb) ~ 1, jura, jura.grid, model = m) : erro na avaliação do argumento 'locations' na seleção do método para a função 'krige': 'Error: object 'jura' not found ' Alguém saberia me explicar qual é o problema? Desde já agradeço, Hélder Gramacho dos Santos Engenheiro Agrônomo Mestrando Ciênc. Geod. Tec. Geoinformação-UFPE *agrohelder(a)gmail.com <agrohelder(a)hotmail.com>*

2 4

Gráficos de barra com duas variáveis sendo comparada
by Rodrigo Muniz 03 Nov '13

03 Nov '13

Caros Membros da Lista, Tenho o seguinte grupo de dados: 0 - 10 10 - 20 Mata Nativa 10.046667 7.600000 Citrus 37 12.363333 23.126667 Citrus 37 15.326667 27.580000 Citrus 24 4.380000 8.856667 Citrus 18 10.430000 21.576667 Citrus 31 9.196667 18.036667 Pasto 37 6.690000 12.533333 Eu gostaria de plotar esses dados em um único gráfico de barras, sendo que as barra correspondente a "Mata nativa" nas camadas 0-10 e 10-20 ficassem lado a lado, e o mesmo para as outras classe. Alguém poderia me dar uma ajuda? Desde já agradeço. -- *Rodrigo A. Muniz* Eng. Agrônomo. Ms Produção Vegetal (UENF) Doutorando em Engenharia de Sistemas Agrícolas (ESALQ/USP) Cel (19) 8300-4333 (Pessoal) Tel (19) 3375-1196 (Residêncial) Tel (19) 3447-8537 (Trabalho)

3 3

Diferença entre séries temporais
by Diogo Caribe de Sousa 03 Nov '13

03 Nov '13

Caros, Vocês podem me dar uma dica de como testo se existe diferença entre duas séries temporais. Obrigado

3 2

Codificação de arquivo .nnw
by Andre Oliveira 02 Nov '13

02 Nov '13

Boa noite, ha 3 anos contrui uma apostila usando o Latex + Sweave .rnw e foi criado no Windows com TinnR. Resolvi atualisar a apostila e compilar o arquivo .rnw, porém, percebi que meu arquivo .rnw está com toda ascentuação comprometida. Onde tem ascentos, til tudo com ponto de interrogação. Recorro a vocês do grupo para tentar resolver este problema. Como mudei para Linux cheguei abrir este arquivo no Linux com Rstudio. Alguma dica? obrigado André Oliveira Souza

1 0

package ‘divisors’ is not available (for R version 3.0.2)
by Andre Oliveira 02 Nov '13

02 Nov '13

Pessoas, preciso da ajuda de vocês estou tentado instalar estes pacotes e obtenho a aviso. install.packages("divisors",dep=T) package ‘divisors’ is not available (for R version 3.0.2) Minha versão consta a informação. R version 3.0.2 (2013-09-25) -- "Frisbee Sailing" Copyright (C) 2013 The R Foundation for Statistical Computing Platform: i686-pc-linux-gnu (32-bit) R é um software livre e vem sem GARANTIA ALGUMA. Você pode redistribuí-lo sob certas circunstâncias. Digite 'license()' ou 'licence()' para detalhes de distribuição. R é um projeto colaborativo com muitos contribuidores. Digite 'contributors()' para obter mais informações e 'citation()' para saber como citar o R ou pacotes do R em publicações. Digite 'demo()' para demonstrações, 'help()' para o sistema on-line de ajuda, ou 'help.start()' para abrir o sistema de ajuda em HTML no seu navegador. Digite 'q()' para sair do R. O que pode estar ocorrendo? obrigado André Oliveira Souza

4 8

[OFF-TOPIC] Arquivos para Tese em Latex UFLA 2013/1
by walmes . 01 Nov '13

01 Nov '13

Saudações, Venho por meio deste divulgar que está disponível em https://github.com/walmes/Tese arquivos para fazer Dissertação e Tese com LaTex nos formatos da UFLA vigentes em 2013/1. Aos que não são acadêmicos da UFLA vale a pena verificar como usar o abntcite para fazer referências bibliográficas segundo as normas da ABNT. Para iniciantes em LaTex, o projeto da Tese inclui um longo preâmbulo para formatação do documento nas normas, além disso tem exemplos de como fazer anexos (inclusive com código R) e rotacionar tabelas que não são tarefas cotidianas, então vale ter como consulta. Outra coisa interessante é a partição do código fonte e a organização em diretórios. À 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 skype: walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================

2 1

Re: [R-br] erro na validação cruzada com pacote geoR
by Hélio Gallo Rocha 31 Out '13

31 Out '13

Caro Eder, Já havia feito a retirada de tendencia com o atributo altura da planta, que está ligado com a produção. O variograma gerado ficou diferente se usar polinomial de 1st, mas não sei como validar isso depois... Abraço Em 31 de outubro de 2013 13:15, Éder Comunello <comunello.eder(a)gmail.com>escreveu: > Boa tarde, Helio! > > Vi seu email, mas vou precisar reservar um tempo pra tratar dele. Você > chegou a usar as glebas como covariável na modelagem? Para isso cada ponto > deveria ter uma var a mais identificando a gleba origem. > > Atte., > > > > Éder Comunello <c <comunello.eder(a)gmail.com>omunello.eder(a)gmail.com> > Dourados, MS - [22 16.5'S, 54 49'W] > > > Em 29 de outubro de 2013 22:25, Hélio Gallo Rocha < > heliogallorocha(a)gmail.com> escreveu: > >> Caro Eder, >> >> Só hoje vi seu post com CRM anexo. >> >> Muito obrigado pela revisão do CRM. >> >> Como também está a procura de uma maior entendimento da geoestatística, >> vamos seguindo. >> >> >> Vou numerar para separar as observações: >> >> 1. Pra entender melhor a área, realmente é café, o recorte esta meio >> estranho, mas na realidade são três talhões de 0.5 ha cada. >> Cada um tem uma característica bem diferenciada com relação a >> produção:alta, baixa e média. >> >> Segundo li, numa situação em que a média é flutuante, como parece ser o >> caso, teríamos de usar krigagem ordinária. >> >> >> A maioria das observações tinha notado, a linha: >> plot(dados, lowess=T, trend='1st') >> achei interessante, mostrando o resultado da retirada de tendencia. >> >> 2. no final : >> >> ### Lembra do resíduo que separei anteriormente? >> ### Se fosse usar... >> va.r = variog(dados.res, uvec=uvec, pairs.min=pairs, max.dist=hmax, >> estimator="classical") >> va.r.env <- variog.mc.env(dados.res, obj.var = va.r) >> plot(va.r, xlab="Dist. (m)", ylim=c(0,summary(va.r$v)[6]*1.1), >> env=va.r.env, main="Residual") >> >> Foi possível gerar o semivariograma sem tendencia, com os resíduos, >> fazer o ajuste assentimento usando: eyefit(va.r) >> >> Dai fazer a validação, fazendo a validação ou não e gerar o mapa sem >> tendencia. >> >> 3. Quanto ao envelopamento, não sei se estou certo, mas ele simula >> valores máximos e minimos para dependencia, mas como é uma simulação, a >> cada vez que vc. roda, pode ter pontos com dependencia, nem que seje um, e >> as vezes não dá nenhum, mas de qualquer forma indica que o valor do alcance >> é pequeno em relação ao todo e baixa dependencia espacial, como nos meus >> dados. >> >> Esse item do envelopamento vale discussão, principalmente em situações >> que a variancia sobe além do limite superior e desce, isso significa algo >> importante, mas o que? >> >> 4. Quanto aos links do Landim já tinha visto, usando o surfer, também >> vale atenção. >> Fiz o recomendado por ele, mas dobrou os valores do mapa., vc. tentou? >> >> 5. Quanto ao ajuste por verosimilhança colocada pelo Prof. Paulo: >> Na primeira vez que havia rodado os dados achei muito estranho, pois o >> semivariograma subia indefinidamente e o semivariograma teórico por >> verossimilhança, usando likfit, baixou muito o patamar. >> Interessante se usar likfit com trend="1st", o semivariograma teórico >> baixa mais ainda, das dai não consigo krigar >> >> 6. Não sei se posso fazer o ajuste assentimento usando o semivariograma >> com tendencia mas usar os parametros sem tendencia. >> >> 7. Tenho de retirar a tendencia sempre? >> Tem casos que há uma patamar claro, mas se rodar o variog com trend, fica >> efeito pepita puro...e ai? >> >> 8. Outra coisa que não sei se estou entendendo: >> após a retirada da tendencia, fazer o mapa com o residuo, mas os valores >> não tem nada a ver com os reais, onde estou entendendo errado? >> >> Abraço, >> >> Hélio >> >> >> >> >> >> Em 26 de outubro de 2013 11:05, Eder Comunello [via R-br] < >> ml-node+s2285057n4660724h12(a)n4.nabble.com> escreveu: >> >>> Senhores, bom dia! >>> >>> Agradeço novamente ao Prof. Paulo pela atenção dispensada. >>> >>> Hélio, quanto ao script que eu tinha ficado de olhar, vou colocar minhas >>> considerações abaixo. >>> >>> Procurei não mudar muito, mas acabei alterando um pouco. Os comentários >>> estão no corpo do script. >>> >>> Lembrando sempre que também estou aprendendo sobre o assunto, então pode >>> ter muita coisa errada ou mal interpretada aí. Mas acho que pode te ajudar. >>> >>> Outra coisa é que não tenho muito ideia do que os dados são, como foram >>> obtidos e as condições da área de onde vieram, o que pra mim é fundamental >>> pra análise. Pelos indícios nos dados (localização, estrutura), vou >>> arriscar que são dados da produção de café ou citros e que o terreno está >>> numa pendente da direita pra esquerda, com maior produção na parte baixa. O >>> que chamou a atenção é o formato da área que é bem particular e não parece >>> muito 'natural'. >>> >>> Mas como o foco principal é o uso do R/geoR, vou tentar me deter nisso. >>> >>> Use por sua própria conta e risco! :D >>> >>> >>> ### <BEGIN> >>> >>> require(geoR) #;require(tcltk2) >>> setwd("C:/LAB/RGIS/trend"); dir() >>> par.ori <- par(no.readonly = TRUE) >>> >>> url1 <- 'https://dl.dropboxusercontent.com/u/117618178/helio/dados.txt' >>> url2 <- 'https://dl.dropboxusercontent.com/u/117618178/helio/recorte.txt >>> ' >>> >>> ### Como criar uma pasta pública no Dropbox... >>> browseURL('https://www.dropbox.com/help/16/pt_BR') >>> >>> ### 1. Carregando dados >>> ------------------------------------------------------ >>> dados <- read.geodata(url1, coords.col=1:2, data.col=3, sep="", header=T) >>> dados$borders <- read.table(url2, header=T) >>> summary(dados); length(dados$data) >>> ls(dados);ls.str(dados) >>> >>> ### 2. Visualizando dados >>> ---------------------------------------------------- >>> >>> ### Visão geral >>> par(mfrow=c(1,1)); par()$mfrow >>> points(dados) >>> points(dados, pt.div="quart") >>> points(dados, pt.div="quint") >>> ### tem pelos três sub-áreas bem marcadas... >>> >>> ### Verificando tendências com o plot.geodata >>> -------------------------------- >>> ### Antes de partir pros variogramas, use o plot.geodata para >>> pré-avaliar tendências >>> ### Bom para esclarecer opções de trend, pois notei certa confusão >>> (cte,1st,2nd) >>> ### Recomendo ler o item 'trend' na ajuda da função plot.geodata >>> ?plot.geodata >>> >>> plot(dados, lowess=T) >>> plot(dados, lowess=T, trend='cte') ### mesmo que anterior >>> ### linhas de comando anteriores são equivalentes: trend='cte' é >>> 'default' da geoR! >>> ### nos gráf. 2 e 3, plota-se as coord vs. dados. Observar as distorções >>> na linha 'lowess' >>> >>> plot(dados, lowess=T, trend=~coords) >>> plot(dados, lowess=T, trend='1st') >>> plot(dados, lowess=T, trend=~dados$coords[,1]+dados$coords[,2]) >>> ### as três linhas anteriores são equivalentes >>> ### o uso da polinomial de primeira ordem melhora o comportamento do >>> histograma >>> ### mas a linha de lowess continua com comportamento anômalo >>> >>> plot(dados, lowess=T, trend='2nd') >>> plot(dados, lowess=T, >>> trend=~I(coords)+I(coords^2)+I(coords[,1]*coords[,2])) >>> ### as duas linhas anteriores são equivalentes (polinomal de segunda >>> ordem)! >>> ### não aparenta agregar nenhuma melhoria >>> >>> ### Mais sobre modelos de tendência com trend.spatial >>> ------------------------ >>> ### Pra melhorar o entendimento dos modelos, dá pra dar uma olhada nas >>> matrizes >>> ?trend.spatial >>> head(unclass(trend.spatial('cte', dados))) ### constante >>> head(unclass(trend.spatial(~coords, dados))) ### polinomal 1a ordem, >>> usando coords >>> head(unclass(trend.spatial('1st', dados))) ### mesmo que anterior >>> head(unclass(trend.spatial('2nd', dados))) ### polinomal 2a ordem, >>> usando coords >>> head(unclass(trend.spatial(~coords[,1], dados))) ### em função de X >>> head(unclass(trend.spatial(~coords[,2], dados))) ### em função de Y >>> head(unclass(trend.spatial(~coords+I(coords[,1]^2), dados))) ### coords >>> + X**2 >>> head(unclass(trend.spatial(~I(coords)+I(coords^2), dados))) >>> head(unclass(trend.spatial(~I(coords)+I(coords^2)+I(coords[,1]*coords[,2]), >>> dados))) ### mesmo que 2nd >>> head(unclass(trend.spatial(~I(coords)+I(coords[,1]^2), dados))) >>> head(unclass(trend.spatial(~I(coords^2), dados))) >>> head(unclass(trend.spatial(~coords[,1]*coords[,2], dados))) >>> head(unclass(trend.spatial(~I(coords[,1]*coords[,2]), dados))) >>> head(unclass(trend.spatial(~coords[,1]+I(coords[,1]^2)+coords[,2], >>> dados))) >>> >>> ### Separando o drift com base na matriz trend.spatial >>> ----------------------- >>> trend.matrix <- unclass(trend.spatial(trend='1st', geodata=dados)); >>> nrow(trend.matrix) >>> trend.lm <- lm(dados$data ~ trend.matrix) >>> trend.res <- lm(dados$data ~ trend.matrix)$residuals >>> trend.fit <- lm(dados$data ~ trend.matrix)$fitted.values >>> >>> ### 'Empacotando' os resíduos em um >>> geodata----------------------------------- >>> dados.res <- dados >>> names(trend.res) <- NULL; trend.res >>> dados.res$data <- trend.res; dados.res$data >>> str(dados.res) >>> ### Pode usar dados.res pra fazer a análise variográfica no lugar de >>> dados, >>> ### caso tenha problemas com krige.conv()!!! >>> >>> ### Embora o uso da polinomial de 1a ordem tenha representado alguma >>> melhoria, >>> ### eu não fiquei muito entusiamado por conta do 'lowess' apresentado >>> ### Talvez fosse o caso de testar covariáveis. >>> ### Uma fácil de obter e que pode ter relação é a altitude... >>> ### Mas vamos investigar o efeito nos variogramas >>> >>> ### 3. Distancia máxima e parâmetros para semivariograma >>> --------------------- >>> hmax <- summary(dados)[[3]][[2]]*.8 ; hmax ### 80% da máxima distância >>> uvec <- 12; pairs <- 20 ### pairs.min >>> >>> ### 4. Variogramas >>> ----------------------------------------------------------- >>> ### Não precisaria rodar todos, afinal como visto, tá repetindo >>> modelos!!! >>> ### v0:média constante; v1:1st; v2:2nd; v3:cte; v4:coords >>> v0 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, >>> estimator="classical") >>> v1 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, >>> estimator="classical", trend="1st") >>> v2 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, >>> estimator="classical", trend="2nd") >>> v3 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, >>> estimator="classical", trend="cte") >>> v4 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, >>> estimator="classical", trend=~coords) >>> >>> ### envelopes >>> v0.env <- variog.mc.env(dados, obj.var = v0) >>> v1.env <- variog.mc.env(dados, obj.var = v1) >>> v2.env <- variog.mc.env(dados, obj.var = v2) >>> v3.env <- variog.mc.env(dados, obj.var = v3) >>> v4.env <- variog.mc.env(dados, obj.var = v4) >>> >>> ### 5. Visualizando variogramas >>> ---------------------------------------------- >>> par(par.ori) >>> par(mfrow=c(2,3)) >>> >>> points(dados, pt.div="quartiles",xlab="Leste", ylab="Norte") >>> >>> plot(v0, xlab="Dist. (m)", ylim=c(0,summary(v0$v)[6]*1.1), env=v0.env, >>> main="Isotrópico com tendência") >>> >>> text(v0$u,v0$v,round(v0$n,1),pos=1);text(v0$u,v0$v,round(v0$u,1),pos=3);lines(v0$u,v0$v);abline(h=mean(v0$v)) >>> >>> plot(v1, xlab="Dist. (m)", ylim=c(0,summary(v1$v)[6]*1.1), env=v1.env, >>> main="Isotrópico trend 1st") >>> >>> text(v1$u,v1$v,round(v1$n,1),pos=1);text(v1$u,v1$v,round(v1$u,1),pos=3);lines(v1$u,v1$v);abline(h=mean(v1$v)) >>> >>> plot(v2, xlab="Dist. (m)", ylim=c(0,summary(v2$v)[6]*1.1), env=v2.env, >>> main="Isotrópico trend 2nd") >>> >>> text(v2$u,v2$v,round(v2$n,1),pos=1);text(v2$u,v2$v,round(v2$u,1),pos=3);lines(v2$u,v2$v);abline(h=mean(v2$v)) >>> >>> plot(v3, xlab="Dist. (m)", ylim=c(0,summary(v3$v)[6]*1.1), env=v3.env, >>> main="Isotrópico trend cte") >>> >>> text(v3$u,v3$v,round(v3$n,1),pos=1);text(v3$u,v3$v,round(v3$u,1),pos=3);lines(v3$u,v3$v);abline(h=mean(v3$v)) >>> >>> plot(v4, xlab="Dist. (m)", ylim=c(0,summary(v4$v)[6]*1.1), env=v4.env, >>> main="Isotrópico trend coords") >>> >>> text(v4$u,v4$v,round(v4$n,1),pos=1);text(v4$u,v4$v,round(v4$u,1),pos=3);lines(v4$u,v4$v);abline(h=mean(v4$v)) >>> >>> ### Analisando os variogramas >>> ------------------------------------------------ >>> ### Realmente o comportamento do histograma de trend='cte' parece conter >>> uma tendência, >>> ### mas os modelos de 'retirada' usados deixam os dados sem uma >>> estrutura de >>> ### variação. Dificilmente vai conseguir ajustar um bom modelo nesses >>> casos, pois >>> ### o comportamento é muito errático. É possível que O erro que você >>> estava obtendo >>> ### seja decorrente da fragilidade do modelo ajustado. >>> >>> ### Realmente não ficou bom! Mas e aí? O que dá pra fazer? Tenta ssim >>> mesmo? >>> ### Ainda tem que testar a anisotropia, que não tem no seu script >>> original. >>> ?variog4 >>> >>> va <- variog4(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, >>> estimator="classical") >>> par(par.ori) >>> plot(va) >>> >>> ### Pronto! Apareceu anisotropia... >>> ### Tem duas sugestões de leitura abaixo pra dar uma ajuda... >>> browseURL(' >>> http://www.rc.unesp.br/igce/aplicada/DIDATICOS/LANDIM/tkrigagem.pdf') >>> browseURL(' >>> https://www.soils.org/publications/sssaj/abstracts/61/1/SS0610010298') >>> >>> ### Você vai ter que estudar adequadamente a questão da anisotropia em >>> seus dados >>> ### Isso não te isenta do efeito combinado anisotropia + tendência >>> ### Só pra testar as possibilidades, vamos fazer o variograma na direção >>> 135 graus >>> >>> ### Vai ter que entrar a direções, sendo mais fácil usar PI/n >>> ### Tabela das direções (PI/n) para variogramas >>> denom <- c(1:4,6,8,10,12,24,36,45,60); k <- (pi/180) ### Denominadores >>> data.frame(rad=paste0('pi/', denom), rad2=round(pi/denom,4), >>> deg=pi/denom/k) >>> ### 135 é 3*45 ou 3*pi/4 >>> >>> va.135 <- variog(dados, uvec=uvec, pairs.min=pairs, dir=3*pi/4, >>> max.dist=hmax, estimator="classical") >>> plot(va.135) >>> ### Aqui parece ser mais fácil ajustar. Talvez tenha que ajustar o 'lag' >>> ### Avaliar a questão da distância entre amostras... >>> >>> ### Lembra do resíduo que separei anteriormente? >>> ### Se fosse usar... >>> va.r = variog(dados.res, uvec=uvec, pairs.min=pairs, max.dist=hmax, >>> estimator="classical") >>> va.r.env <- variog.mc.env(dados.res, obj.var = va.r) >>> plot(va.r, xlab="Dist. (m)", ylim=c(0,summary(va.r$v)[6]*1.1), >>> env=va.r.env, main="Residual") >>> >>> ### Vou parar por aqui, porque antes de prosseguir vai precisar ver >>> ### adequadamente essa questão da anisotropia. >>> ### No caso do efeito ser só da anisotropia, vai simplificar a krigagem >>> e você >>> ### não deve mais ter o problema que vinha tendo. >>> >>> ### <END> >>> >>> Éder Comunello <[hidden email]<http://user/SendEmail.jtp?type=node&node=4660724&i=0>[hidden >>> email] <http://user/SendEmail.jtp?type=node&node=4660724&i=1>> >>> Dourados, MS - [22 16.5'S, 54 49'W] >>> >>> >>> >>> >>> _______________________________________________ >>> R-br mailing list >>> [hidden email] <http://user/SendEmail.jtp?type=node&node=4660724&i=2> >>> 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. >>> >>> ------------------------------ >>> If you reply to this email, your message will be added to the >>> discussion below: >>> >>> http://r-br.2285057.n4.nabble.com/Re-R-br-erro-na-validacao-cruzada-com-pac… >>> To unsubscribe from R-br, click here<http://r-br.2285057.n4.nabble.com/template/NamlServlet.jtp?macro=unsubscrib…> >>> . >>> NAML<http://r-br.2285057.n4.nabble.com/template/NamlServlet.jtp?macro=macro_view…> >>> >> >> >> >> -- >> Hélio Gallo Rocha >> IFSULDEMINAS - Câmpus Muzambinho >> >> _______________________________________________ >> 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. >> > > -- Hélio Gallo Rocha IFSULDEMINAS - Câmpus Muzambinho

1 0

Legenda em gráfico de regressão
by Edson Lira 31 Out '13

31 Out '13

Caros amigos, tenho a rotina abaixo: require(PASWR) data(vit2005) attach(vit2005) names(vit2005) head(vit2005) Elevator<-as.factor(elevator) contrasts(Elevator) modSimpl<-lm(totalprice~area) summary(modSimpl) plot(totalprice~area) Este banco tem uma variável que é o Elevator (presença ou ausência do elevador), gostaria de plotar o gráfico da regressão totalprice~area, identificando os pontos aonde tem a presença ou ausência do elevador. Como faço? plot(totalprice~area) Edson Lira Estatístico Manaus-Amazonas

3 2

ggplot2
by geovane barbosa 30 Out '13

30 Out '13

OLá pessoal tudo bem, estou usando o seguinte código e esta dando o seguinte erro Error: Aesthetics must either be length one, or the same length as the dataProblems:IDADEA, factor(Cor2) Alguém poderia me dar uma dica aguardo #################################### #################################### Prof. Geovane Carlos Barbosa UCL - Faculdade do Centro Leste ################################### ###################################

3 2

Variável Dummy na regressão
by Edson Lira 29 Out '13

29 Out '13

Caros amigos, estou tentando fazer um estudo de regressão com duas variáveis preditoras. Dados e comandos abaixo: mus<-c(2,6,8,3,2,7,9,8,4,6) tr<-c("e","m","e","m","e","e","m","m","e","e") trh<-c(2.9,3.0,4.8,1.8,2.9,4.9,4.2,4.8,4.4,4.5) tipo<-ifelse(tr=="e",1,0) #codificar a variável tr, 1=eletrico 2=mecânico b1<-data.frame(mus,tipo,trh) Quero fazer a análise usando a variável tipo. São dois modelos: um quando o e=1 e o outro e=0. m1<-lm(trh~mus+tipo=0, data=b1) m2<-lm(trh~mus+tipo=1, data=b1) Não estou conseguindo implementar, alguém poderia dar alguma sugestão? [ ]'s. Edson Lira Estatístico Manaus-Amazonas

4 4

Re: [R-br] erro na validação cruzada com pacote geoR
by Éder Comunello 29 Out '13

29 Out '13

Senhores, bom dia! Agradeço novamente ao Prof. Paulo pela atenção dispensada. Hélio, quanto ao script que eu tinha ficado de olhar, vou colocar minhas considerações abaixo. Procurei não mudar muito, mas acabei alterando um pouco. Os comentários estão no corpo do script. Lembrando sempre que também estou aprendendo sobre o assunto, então pode ter muita coisa errada ou mal interpretada aí. Mas acho que pode te ajudar. Outra coisa é que não tenho muito ideia do que os dados são, como foram obtidos e as condições da área de onde vieram, o que pra mim é fundamental pra análise. Pelos indícios nos dados (localização, estrutura), vou arriscar que são dados da produção de café ou citros e que o terreno está numa pendente da direita pra esquerda, com maior produção na parte baixa. O que chamou a atenção é o formato da área que é bem particular e não parece muito 'natural'. Mas como o foco principal é o uso do R/geoR, vou tentar me deter nisso. Use por sua própria conta e risco! :D ### <BEGIN> require(geoR) #;require(tcltk2) setwd("C:/LAB/RGIS/trend"); dir() par.ori <- par(no.readonly = TRUE) url1 <- 'https://dl.dropboxusercontent.com/u/117618178/helio/dados.txt' url2 <- 'https://dl.dropboxusercontent.com/u/117618178/helio/recorte.txt' ### Como criar uma pasta pública no Dropbox... browseURL('https://www.dropbox.com/help/16/pt_BR') ### 1. Carregando dados ------------------------------------------------------ dados <- read.geodata(url1, coords.col=1:2, data.col=3, sep="", header=T) dados$borders <- read.table(url2, header=T) summary(dados); length(dados$data) ls(dados);ls.str(dados) ### 2. Visualizando dados ---------------------------------------------------- ### Visão geral par(mfrow=c(1,1)); par()$mfrow points(dados) points(dados, pt.div="quart") points(dados, pt.div="quint") ### tem pelos três sub-áreas bem marcadas... ### Verificando tendências com o plot.geodata -------------------------------- ### Antes de partir pros variogramas, use o plot.geodata para pré-avaliar tendências ### Bom para esclarecer opções de trend, pois notei certa confusão (cte,1st,2nd) ### Recomendo ler o item 'trend' na ajuda da função plot.geodata ?plot.geodata plot(dados, lowess=T) plot(dados, lowess=T, trend='cte') ### mesmo que anterior ### linhas de comando anteriores são equivalentes: trend='cte' é 'default' da geoR! ### nos gráf. 2 e 3, plota-se as coord vs. dados. Observar as distorções na linha 'lowess' plot(dados, lowess=T, trend=~coords) plot(dados, lowess=T, trend='1st') plot(dados, lowess=T, trend=~dados$coords[,1]+dados$coords[,2]) ### as três linhas anteriores são equivalentes ### o uso da polinomial de primeira ordem melhora o comportamento do histograma ### mas a linha de lowess continua com comportamento anômalo plot(dados, lowess=T, trend='2nd') plot(dados, lowess=T, trend=~I(coords)+I(coords^2)+I(coords[,1]*coords[,2])) ### as duas linhas anteriores são equivalentes (polinomal de segunda ordem)! ### não aparenta agregar nenhuma melhoria ### Mais sobre modelos de tendência com trend.spatial ------------------------ ### Pra melhorar o entendimento dos modelos, dá pra dar uma olhada nas matrizes ?trend.spatial head(unclass(trend.spatial('cte', dados))) ### constante head(unclass(trend.spatial(~coords, dados))) ### polinomal 1a ordem, usando coords head(unclass(trend.spatial('1st', dados))) ### mesmo que anterior head(unclass(trend.spatial('2nd', dados))) ### polinomal 2a ordem, usando coords head(unclass(trend.spatial(~coords[,1], dados))) ### em função de X head(unclass(trend.spatial(~coords[,2], dados))) ### em função de Y head(unclass(trend.spatial(~coords+I(coords[,1]^2), dados))) ### coords + X**2 head(unclass(trend.spatial(~I(coords)+I(coords^2), dados))) head(unclass(trend.spatial(~I(coords)+I(coords^2)+I(coords[,1]*coords[,2]), dados))) ### mesmo que 2nd head(unclass(trend.spatial(~I(coords)+I(coords[,1]^2), dados))) head(unclass(trend.spatial(~I(coords^2), dados))) head(unclass(trend.spatial(~coords[,1]*coords[,2], dados))) head(unclass(trend.spatial(~I(coords[,1]*coords[,2]), dados))) head(unclass(trend.spatial(~coords[,1]+I(coords[,1]^2)+coords[,2], dados))) ### Separando o drift com base na matriz trend.spatial ----------------------- trend.matrix <- unclass(trend.spatial(trend='1st', geodata=dados)); nrow(trend.matrix) trend.lm <- lm(dados$data ~ trend.matrix) trend.res <- lm(dados$data ~ trend.matrix)$residuals trend.fit <- lm(dados$data ~ trend.matrix)$fitted.values ### 'Empacotando' os resíduos em um geodata----------------------------------- dados.res <- dados names(trend.res) <- NULL; trend.res dados.res$data <- trend.res; dados.res$data str(dados.res) ### Pode usar dados.res pra fazer a análise variográfica no lugar de dados, ### caso tenha problemas com krige.conv()!!! ### Embora o uso da polinomial de 1a ordem tenha representado alguma melhoria, ### eu não fiquei muito entusiamado por conta do 'lowess' apresentado ### Talvez fosse o caso de testar covariáveis. ### Uma fácil de obter e que pode ter relação é a altitude... ### Mas vamos investigar o efeito nos variogramas ### 3. Distancia máxima e parâmetros para semivariograma --------------------- hmax <- summary(dados)[[3]][[2]]*.8 ; hmax ### 80% da máxima distância uvec <- 12; pairs <- 20 ### pairs.min ### 4. Variogramas ----------------------------------------------------------- ### Não precisaria rodar todos, afinal como visto, tá repetindo modelos!!! ### v0:média constante; v1:1st; v2:2nd; v3:cte; v4:coords v0 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical") v1 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical", trend="1st") v2 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical", trend="2nd") v3 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical", trend="cte") v4 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical", trend=~coords) ### envelopes v0.env <- variog.mc.env(dados, obj.var = v0) v1.env <- variog.mc.env(dados, obj.var = v1) v2.env <- variog.mc.env(dados, obj.var = v2) v3.env <- variog.mc.env(dados, obj.var = v3) v4.env <- variog.mc.env(dados, obj.var = v4) ### 5. Visualizando variogramas ---------------------------------------------- par(par.ori) par(mfrow=c(2,3)) points(dados, pt.div="quartiles",xlab="Leste", ylab="Norte") plot(v0, xlab="Dist. (m)", ylim=c(0,summary(v0$v)[6]*1.1), env=v0.env, main="Isotrópico com tendência") text(v0$u,v0$v,round(v0$n,1),pos=1);text(v0$u,v0$v,round(v0$u,1),pos=3);lines(v0$u,v0$v);abline(h=mean(v0$v)) plot(v1, xlab="Dist. (m)", ylim=c(0,summary(v1$v)[6]*1.1), env=v1.env, main="Isotrópico trend 1st") text(v1$u,v1$v,round(v1$n,1),pos=1);text(v1$u,v1$v,round(v1$u,1),pos=3);lines(v1$u,v1$v);abline(h=mean(v1$v)) plot(v2, xlab="Dist. (m)", ylim=c(0,summary(v2$v)[6]*1.1), env=v2.env, main="Isotrópico trend 2nd") text(v2$u,v2$v,round(v2$n,1),pos=1);text(v2$u,v2$v,round(v2$u,1),pos=3);lines(v2$u,v2$v);abline(h=mean(v2$v)) plot(v3, xlab="Dist. (m)", ylim=c(0,summary(v3$v)[6]*1.1), env=v3.env, main="Isotrópico trend cte") text(v3$u,v3$v,round(v3$n,1),pos=1);text(v3$u,v3$v,round(v3$u,1),pos=3);lines(v3$u,v3$v);abline(h=mean(v3$v)) plot(v4, xlab="Dist. (m)", ylim=c(0,summary(v4$v)[6]*1.1), env=v4.env, main="Isotrópico trend coords") text(v4$u,v4$v,round(v4$n,1),pos=1);text(v4$u,v4$v,round(v4$u,1),pos=3);lines(v4$u,v4$v);abline(h=mean(v4$v)) ### Analisando os variogramas ------------------------------------------------ ### Realmente o comportamento do histograma de trend='cte' parece conter uma tendência, ### mas os modelos de 'retirada' usados deixam os dados sem uma estrutura de ### variação. Dificilmente vai conseguir ajustar um bom modelo nesses casos, pois ### o comportamento é muito errático. É possível que O erro que você estava obtendo ### seja decorrente da fragilidade do modelo ajustado. ### Realmente não ficou bom! Mas e aí? O que dá pra fazer? Tenta ssim mesmo? ### Ainda tem que testar a anisotropia, que não tem no seu script original. ?variog4 va <- variog4(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical") par(par.ori) plot(va) ### Pronto! Apareceu anisotropia... ### Tem duas sugestões de leitura abaixo pra dar uma ajuda... browseURL(' http://www.rc.unesp.br/igce/aplicada/DIDATICOS/LANDIM/tkrigagem.pdf') browseURL(' https://www.soils.org/publications/sssaj/abstracts/61/1/SS0610010298') ### Você vai ter que estudar adequadamente a questão da anisotropia em seus dados ### Isso não te isenta do efeito combinado anisotropia + tendência ### Só pra testar as possibilidades, vamos fazer o variograma na direção 135 graus ### Vai ter que entrar a direções, sendo mais fácil usar PI/n ### Tabela das direções (PI/n) para variogramas denom <- c(1:4,6,8,10,12,24,36,45,60); k <- (pi/180) ### Denominadores data.frame(rad=paste0('pi/', denom), rad2=round(pi/denom,4), deg=pi/denom/k) ### 135 é 3*45 ou 3*pi/4 va.135 <- variog(dados, uvec=uvec, pairs.min=pairs, dir=3*pi/4, max.dist=hmax, estimator="classical") plot(va.135) ### Aqui parece ser mais fácil ajustar. Talvez tenha que ajustar o 'lag' ### Avaliar a questão da distância entre amostras... ### Lembra do resíduo que separei anteriormente? ### Se fosse usar... va.r = variog(dados.res, uvec=uvec, pairs.min=pairs, max.dist=hmax, estimator="classical") va.r.env <- variog.mc.env(dados.res, obj.var = va.r) plot(va.r, xlab="Dist. (m)", ylim=c(0,summary(va.r$v)[6]*1.1), env=va.r.env, main="Residual") ### Vou parar por aqui, porque antes de prosseguir vai precisar ver ### adequadamente essa questão da anisotropia. ### No caso do efeito ser só da anisotropia, vai simplificar a krigagem e você ### não deve mais ter o problema que vinha tendo. ### <END> Éder Comunello <c <comunello.eder(a)gmail.com>omunello.eder(a)gmail.com> Dourados, MS - [22 16.5'S, 54 49'W]

2 1

Recusão no R, problemas ao implementar o algoritimo insertionsort no R
by Augusto Ribas 29 Out '13

29 Out '13

Ola pessoa, eu estava tentando implementar o algorítimo insertion sort no R de forma recursiva, mas não estou conseguindo e também não consigo identificar onde estou errado. #Veja so. Eu tenho um vetor qualquer: vetor<-sample(100) vetor #Iterativamente ele funciona beleza, eu fiz assim: insertionsort<-function(vetor){ n<-length(vetor) for(i in 2:n) { aux=vetor[i] j=i-1 while(j>=1 && vetor[j]>aux) { vetor[j+1]<-vetor[j] j=j-1 } vetor[j+1]=aux } return(vetor) } insertionsort(vetor) #Mas quando tento usar recursão, não consigo fazer funcionar, ele não organiza insertionsortR<-function(vetor,n) { if(n>1) { insertionsortR(vetor,n-1) aux=vetor[n] i=n while(vetor[i-1]>aux && i>1) { vetor[i]<-vetor[i-1] i<-i-1 } vetor[i]<-aux } return(vetor) } insertionsortR(vetor,length(vetor)) #Usando o rcpp ele funcionou dessa forma aqui: cppFunction(" NumericVector insertionsortRC(NumericVector vetor, int n) { double aux; int i; if(n>1) { insertionsortRC(vetor,n-1); aux=vetor[n-1]; i=n-1; while(vetor[i-1]>aux && i>=0 ) { vetor[i]=vetor[i-1]; i--; } vetor[i]=aux; } return vetor; } ") insertionsortRC(vetor,length(vetor)) Alguém sabe dizer onde estou errado la no insertionsort no R. Eu imagino se estou imaginado errado o funcionamento dos vetores no R, se ele ta refazendo um novo vetor a cada chamada e no final nada da certo do jeito que estou fazendo. Recursão funciona beleza no R, por exemplo usando fatorial: fatorial<-function(n) { if(n==1) { return(1) } else { return(n*fatorial(n-1)) } } fatorial(5) Mas no caso do vetor, não entendi o funcionamento ali dos vetores, se alguém puder dar uma luz :) -- Grato Augusto C. A. Ribas Site Pessoal: http://recologia.com.br/ <http://augustoribas.heliohost.org> Github: https://github.com/Squiercg Lattes: http://lattes.cnpq.br/7355685961127056

2 2

Selecionar o primeiro numero negativo depois do último positivo.
by Alisson Lucrecio 29 Out '13

29 Out '13

Caro Colegas da lista r-br, Boa tarde. Eu estou precisando encontrar o primeiro numero negativo depois do último numero positivo de uma sequencia c(10, 12, 11, 5, 4, -5, -10, -11, -15, 4, 5, 6, 8, 10, 8, 5, 4, -6, -2, 4, 8). Assim o numero selecionado seriam -5 e -6. Algum de vocês poderia me ajudar com uma função que tivesse tal finalidade? Obrigado. Alisson Lucrécio da Costa

3 3

Não consigo fazer uma ANOVA com arquivo com 1e+06 registros.
by Luiz Roberto Martins Pinto 29 Out '13

29 Out '13

Caros companheiros da R-BR. Não consigo fazer uma ANOVA com arquivo com 1e+06 registros. Então... preciso de ajuda!!! Dados: http://www1.datafilehost.com/d/c0d31775 Meu pc R for windows 2.15.1(x64) 8 Gb de Memo load('RCBD_Data.Rdata') # Arquivo com 1e+06 registros m=aov(y~Treat+block,data=Data) summary(m) Mensagem de erro: Erro: não é possível alocar vetor de tamanho 14.9 Gb Além disso: Mensagens de aviso perdidas: 1: In model.matrix.default(mt, mf, contrasts) : Reached total allocation of 8086Mb: see help(memory.size) 2: In model.matrix.default(mt, mf, contrasts) : Reached total allocation of 8086Mb: see help(memory.size) 3: In model.matrix.default(mt, mf, contrasts) : Reached total allocation of 8086Mb: see help(memory.size) 4: In model.matrix.default(mt, mf, contrasts) : Reached total allocation of 8086Mb: see help(memory.size) 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

6 24

Problemas com a atualização e instalação de pacotes no R no UBUNTU 13.04
by Luciano F. Sgarbi 29 Out '13

29 Out '13

Prezados, Sou novo no uso do OS linux e estou com problemas na atualização do R (versão 3.0.2) e na instalação de qualquer pacote no R. Meu sistema operacional é o UBUNTU 13.04, 32-bit. Se alguém puder dar uma força serei grato. Abraços, Luciano -- *Luciano F. Sgarbi * * Currículo Lattes <http://lattes.cnpq.br/9788198135883354> **Dep. Ecologia, Instituto de Ciências Biológicas I** * *Universidade Federal de Goiás, campus II, Goiânia-GO, Brasil*

2 5

Polinômios ortogonais
by Fernando Souza 28 Out '13

28 Out '13

Caros amigos, estou realizando uma regressão linear com contrastes polinomiais par obter o grau do polinômio que melhor se ajusta aos dados. Realizei o modelo no R e fiz também manualmente utilizando a tabela de divisores, multiplicadores para comparações ortogonais. Acontece que embora a análise de variância indique ser sigignificativo o modelo de 2ª ordem quando obtenho o modelo pelo método summary a regressão calculada manualmente difere daquela apresentada pelo R alguem pode me explicar o que está acontecendo. Segue o CMR: dados<-structure(list(Bloco = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L), .Label = c("1", "2", "3", "4", "5", "6"), class = "factor"), Espacamento = c(1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5, 1, 2, 3, 4, 5), Producao = c(3.02, 2.795, 2.966, 2.552, 2.822, 3.334, 3.1, 2.651, 2.687, 2.543, 3.065, 2.741, 2.624, 2.84, 2.597, 3.109, 2.85, 2.759, 2.903, 2.57, 3.181, 2.759, 2.759, 2.525, 2.66, 3.244, 2.723, 2.507, 2.417, 3.002)), .Names = c("Bloco", "Espacamento", "Producao"), row.names = c(NA, -30L), class = "data.frame") dados$Bloco<-factor(dados$Bloco) dados$Espacamento<-factor(dados$Espacamento) modelo3<-lm(Producao~Bloco+ordered(Espacamento), data=dados) O modelo obtido manualmente foi: Y =4,721 - 0,0458X + 0,00025X² atenciosamente -- ===================================================== Fernando Souza Zootecnista, DSc. Produção Animal cel: (+55) 82 8113-8781 e-mail:nandodesouza@gmail.com ==================================================

3 3

[OFF] Mensagem duplicada!
by Andre Oliveira 27 Out '13

27 Out '13

Pessoal a lista se encontra com algum problema? Minhas duas ultimas postaguens foi enviada duas vezes na mesma mensagem sendo que enviei somente uma vez! André Oliveira Souza

1 0

Instalar pacotes!
by Andre Oliveira 27 Out '13

27 Out '13

Pessoas, preciso da ajuda de vocês estou tentado instalar estes pacotes e obtenho a aviso. package ‘divisors’ is not available (for R version 3.0.2) install.packages("descr"), dep=T) install.packages("itm", dep=T) install.packages("divisors",dep=T) Se eles não estão disponíveis para o R, minha versão consta a informação. andre@andre-Aspire-5920:~$ R -version ATEN??O: op??o desconhecida '-version' R version 3.0.2 (2013-09-25) -- "Frisbee Sailing" Copyright (C) 2013 The R Foundation for Statistical Computing Platform: i686-pc-linux-gnu (32-bit) R é um software livre e vem sem GARANTIA ALGUMA. Você pode redistribuí-lo sob certas circunstâncias. Digite 'license()' ou 'licence()' para detalhes de distribuição. R é um projeto colaborativo com muitos contribuidores. Digite 'contributors()' para obter mais informações e 'citation()' para saber como citar o R ou pacotes do R em publicações. Digite 'demo()' para demonstrações, 'help()' para o sistema on-line de ajuda, ou 'help.start()' para abrir o sistema de ajuda em HTML no seu navegador. Digite 'q()' para sair do R. [Área de trabalho anterior carregada] O que devo proceder para resolver? obrigado André Oliveira Souza

1 0

Ler banco de dados .xls
by Andre Oliveira 27 Out '13

27 Out '13

Boa tarde, eu sempre usei no linux os pacotes. 1) gdata 2) RODBC Existe algum outro pacote mais moderno e eficiente no momento? André Oliveira Souza

3 2

média em um "replicate"
by Maurício Lordêlo 26 Out '13

26 Out '13

Caros, Devido a um problema com a função "bild" (pertencente ao pacote de mesmo nome), tive que modificar um script de simulação e não estou conseguindo finalizá-lo. No CRM abaixo, a função "geral" ajusta o modelo e extrai os coeficientes de interesse. Em seguida, coloco um "for" com o número de simulações que desejo (coloquei 5 mas pretendo 2000). Ao usar o "replicate" obtenho os resultados. Preciso obter uma média por coluna destes resultados. Pelo que percebi, o objeto que denominei de "repeticoes" armazena um objeto em cada uma de suas linhas. Alguém saberia uma forma de solucionar isso, ou seja, encontrar uma média por coluna do objeto "repeticoes"? Caso não seja esta a opção adequada, aceito sugestões. Agradeço desde já. library(bild) nt <- 8; ntr <- 3 ; nr <- 5 n <- ntr*nr*nt; nid <- ntr * nr id <- rep(1:nid, each = nt) tempo <- rep(1:nt, times = nid) trat <- rep(factor(LETTERS[0:ntr]), each = nr*nt) dados=data.frame(id,tempo,trat) geral <- function() { z <- numeric(n) for(i in 1:n) z[i] <- rbinom(1,1,0.5) dat <- data.frame(dados, y=z) (bild(y ~ trat+tempo, data = dat, time="tempo",start=NULL,dependence="MC1R"))@coefficients[1:3] } ns=5 ### número de simulações for(i in 1:ns) { repeticoes <- replicate(1, geral(), simplify = FALSE) print(repeticoes) }

4 6

helio rocha shared "geoR" with you
by helio rocha via Dropbox 25 Out '13

25 Out '13

Helio used Dropbox to share some files with you! Click here to view geoR: https://www.dropbox.com/l/O4tRyduR1R7iTTAZoIeYLa?text=1

2 1

divergência dos resultados na máquina de 32bits e 64bits (ode e optim)
by Adriele Giaretta Biase 25 Out '13

25 Out '13

Olá pessoal, estou trabalhando com resoluções de Equações Diferenciais Ordinárias no R, usando a função ode do pacote (deSolve), método “lsoda” e também estou usando otimização do sistema de EDO usando a função optim do pacote (stats), método “L-BFGS-B”. Quando executo a mesma programação e com o mesmo banco de dados, a predição dos parâmetros apresentam diferenças entre o R de 32 bits e o R de 64 bits. As versões do R são as mesmas (3.01), portanto, os pacotes são os mesmos. Os resultados podem diferir mesmo usando a mesma configuração das máquina, exceto por usar 32 bits e 64 bits? Lembrando, a versão, os pacotes, a configuração da máquina, a programação e os dados são os mesmos. O fato de usar uma máquina de 32 bits e 64 bits interfere nos resultados por estar usando esses algoritmos de aproximações numéricas? -- Adriele Giaretta Biase. Mestre em Estatística e Experimentação Agropecuária - UFLA. Doutoranda em Estatística e Experimentação Agronômica - ESALQ/ USP Contato: (19) 8861-0619.

4 3

Re: [R-br] erro na validação cruzada com pacote geoR
by Éder Comunello 24 Out '13

24 Out '13

Senhores, boa noite! Elias Krainski e Prof. Paulo, agradeço pelos esclarecimentos prestados! Ainda preciso dedicar um tempo para estudar adequadamente o que foi colocado, mas por ora gostaria apenas de esclarecer se a 'retirada' do efeito da tendência (por ajuste polinomial) é uma estratégia válida ou há alguma ressalva a ser feita. Agradeço se os senhores, ou outro colega da lista, puderem ponderar a respeito. Hélio, vi que você disponibilizou o script e pretendo testá-lo o mais breve possível. Espero que você ainda tenha algum tempo pra resolver! :D Sem mais, Éder Comunello <c <comunello.eder(a)gmail.com>omunello.eder(a)gmail.com> Dourados, MS - [22 16.5'S, 54 49'W]

3 4

Composite Reliability
by Felipe Buchbinder 24 Out '13

24 Out '13

Prezados, Alguém conhece algum pacote e função para calcular *Composite Reliability*no R? Grato, Felipe

2 1

Rcitrus
by ASANTOS 24 Out '13

24 Out '13

Boa tarde pessoal, Gostaria de saber do Elias ou do Paulo ou de algum membro que trabalha com o pacote Rcitrus, qual o valor utilizado como parâmetro que é comparado com o índice calculado (razão da variância observada pela variância teórica do modelo binomial) para decidir se a distribuição é aleatória, agregada ou uniforme no índice de dispersão binomial da função disp.quadrats(), apesar de a saída da função já apresentar a distribuição espacial apresentadas pelos dados avaliados. 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

Dúvida pedigreemm vs lmer
by Jane Freitas 23 Out '13

23 Out '13

Mauro, envio outro anexo. Testei aqui e rodou legal. Segue a rotina abaixo. Agora, você ou alguém da lista sabe dizer por que os BLUPs do modelo sem pedigree (Valor genético = aditivo+ dominância + epistático) foi menor que os BLUPs dom modelo com o pedigree (Valor genético aditivo)??? > dadped<-read.table("genealogia2.csv", head=T, sep=";", dec=",") > head(dadped) Num Id Mon Dad 1 1 GPA518 <NA> <NA> 2 2 GX458 <NA> <NA> 3 3 GH3156 GPA518 <NA> 4 4 GS214 GX458 <NA> 5 5 RA1251 GS214 GH3156 6 6 AB1542 RA1251 GPA518 > require(lme4) > require(pedigreemm) > pedCana <- pedigree(sire = as.integer(dadped$Dad), + dam = as.integer(dadped$Mon), + label= as.character(1:9)) > fac <- relfactor(pedCana) > MpedCana <- crossprod(fac) > dat2<-read.table("dados_sem_ascestrais.csv", head=T, sep=";", dec=",", na.strings="NA") > head(dat2) geno id bloco corte local Prod 1 RA1251 5 1 1 1 117.65 2 RA1251 5 2 1 1 125.35 3 RA1251 5 3 1 1 135.02 4 AB1542 6 1 1 1 101.03 5 AB1542 6 2 1 1 108.59 6 AB1542 6 3 1 1 108.45 > attach(dat2) The following object(s) are masked from 'dat2 (position 3)': bloco, corte, geno, id, local, Prod > Bloco1 <- factor(bloco) > Corte1 <- factor(corte) > Local1 <- factor(local) > fm2 <- pedigreemm(Prod ~ (1|id) + Local1/Corte1 + (1|Local1/Bloco1), data = dat2, pedigree = list(id = pedCana)) > BLUP.a <- ranef(fm2)$id; BLUP.a # Valor genético aditivo (Intercept) 5 0.4146977 6 -1.1289438 7 1.2858838 8 0.1372509 9 0.2615716 > fm3 <-lmer(Prod ~ (1|id) + Local1/Corte1 + (1|Local1/Bloco1), data = dat2) > BLUP.g <- ranef(fm3)$id; BLUP.g # Valor genético (Intercept) 5 0.04624036 6 -0.51780417 7 0.30066464 8 0.06359092 9 0.10730825 > data.frame(BLUP.g, BLUP.a) X.Intercept. X.Intercept..1 5 0.04624036 0.4146977 6 -0.51780417 -1.1289438 7 0.30066464 1.2858838 8 0.06359092 0.1372509 9 0.10730825 0.2615716 > -- ================================================= Edjane Gonçalves de Freitas Engenheira Agrônoma - UFAL Doutoranda - Lab. Genética e Estatística Departamento de Genética e Melhoramento de Plantas - ESALQ/USP Cel: (19) 8102-7790 ==================================================

1 0

erro na validação cruzada com pacote geoR
by Hélio Gallo Rocha 22 Out '13

22 Out '13

Caros da lista. Estou uma situação bastante complicada. Estou no fim do mestrado, onde utilizei o R todo o tempo. Mas estou numa situação que talvez tenha de deixar de usar o R e passar para um "pago". Já postei sobre este problema e vi outros com o mesmo. Após a retirada de tendencia do semivariograma e fazer a validação cruzada para escolha do melhor modelo o R dá esse erro: xErro em solve.default(ttivtt, crossprod(ivtt, as.vector(data))) : sistema é computacionalmente singular: condição recíproca número = 1.18641e-024 Alguém poderia me dar uma luz? Não tem solução? é isso? grato a todos -- Hélio Gallo Rocha IFSULDEMINAS - Câmpus Muzambinho

6 15

(sem assunto)
by Jane Freitas 22 Out '13

22 Out '13

Bom dia!! Pessoal, ajustei um modelo utilizando a função pedigreemm (incorporando a informação do pedigree) e lmer (do pacote, lme4, sem o pedigree). Notei que os BLUps obtidos (via ranef) com a pedigreemm são maiores do que os obtidos com a lmer. A dúvida é que acredito que quando usamos a informação do pedigree, predizemos o *Valor Genético Aditivo* (*Breeding Values*). E quando não utilizamos a informação do pedigree, predizemos o *Valor Genético Total* (*incluindo os efeitos aditivos, dominância e epistático*), por isso esperava maiores valores dos BLUPs. Será que alguém pode me esclarecer?? Segue a rotina e dados em anexo: #----------------------------------------------------------------------------- # EXEMPLO: # Variável analisada: Prod # Número de indivíduos no pedigree: 9 # Número de indivíduos com informação fenotipica: 5 # Número de bloco: 3 # Número de corte: 2 # Número de local: 3 #------------------------------------------------------------------------------ rm(list = ls()) *1 - Com o pedigree* *# reading pedigree file* dadped<-read.table("genealogia2.csv", head=T, sep=";", dec=",") head(dadped) require(lme4) require(pedigreemm) *# Constructing pedigree* pedCana <- pedigree(sire = as.integer(dadped$Dad), dam = as.integer(dadped$Mon), label= as.character(1:9)) fac <- relfactor(pedCana) MpedCana <- crossprod(fac) *# reading phenotype file* dat2<-read.table("dados_sem_ascestrais.csv", head=T, sep=";", dec=",", na.strings="NA") head(dat2) attach(dat2) Bloco1 <- factor(bloco) Corte1 <- factor(corte) Local1 <- factor(local) *# fitting model with pedigree* fm2 <- pedigreemm(Prod ~ (1|id) + Local1/Corte1 + (1|Local1/Bloco1), data = dat2, pedigree = list(id = pedCana)) BLUP.a <- ranef(fm2)$id; BLUP.a * # BLUPs dos efeitos genéticos aditivo* *2 - Sem o pedigree* * # fitting model without pedigree* fm3 <-lmer(Prod ~ (1|id) + Local1/Corte1 + (1|Local1/Bloco1), data = dat2) BLUP.g <- ranef(fm3)$id; BLUP.g *# BLUPs dos efeitos genéticos (aditivo + dominância + epistáticos) * data.frame(BLUP.g, BLUP.a) *# BLUPs dos efeitos genéticos e genéticos aditivo* *Resultados:* ** *> fm2* Linear mixed model fit by REML Formula: Prod ~ (1 | id) + Local1/Corte1 + (1 | Local1/Bloco1) Data: dat2 AIC BIC logLik deviance REMLdev 764.3 789.3 -372.1 777.8 744.3 Random effects: Groups Name Variance Std.Dev. Bloco1:Local1 (Intercept) 1.0992e-13 3.3155e-07 id (Intercept) 4.8287e+00 2.1974e+00 Local1 (Intercept) 3.0252e+01 5.5002e+00 Residual 3.3712e+02 1.8361e+01 Number of obs: 90, groups: Bloco1:Local1, 9; id, 5; Local1, 3 *> fm3* Linear mixed model fit by REML Formula: Prod ~ (1 | id) + Local1/Corte1 + (1 | Local1/Bloco1) Data: dat2 AIC BIC logLik deviance REMLdev 764.4 789.4 -372.2 777.8 744.4 Random effects: Groups Name Variance Std.Dev. Bloco1:Local1 (Intercept) 1.1257e-17 3.3551e-09 id (Intercept) 1.3792e+00 1.1744e+00 Local1 (Intercept) 2.9551e+01 5.4361e+00 Residual 3.3926e+02 1.8419e+01 Number of obs: 90, groups: Bloco1:Local1, 9; id, 5; Local1, 3 *> data.frame(BLUP.g, BLUP.a)* X.Intercept. X.Intercept..1 5 0.04624036 0.4146977 6 -0.51780417 -1.1289438 7 0.30066464 1.2858838 8 0.06359092 0.1372509 9 0.10730825 0.2615716 -- ================================================= Edjane Gonçalves de Freitas Engenheira Agrônoma - UFAL Doutoranda - Lab. Genética e Estatística Departamento de Genética e Melhoramento de Plantas - ESALQ/USP Cel: (19) 8102-7790 ==================================================

2 1

Rcpp cppFunction
by Carlos Trucios Maza 22 Out '13

22 Out '13

Pessoal boa tarde, Estou tentando reproducir um código do livro sobre Rcpp e não estou conseguindo. Alguém tem alguma ideia de porque não funciona? (Trabalho com OS X 10.8.4). code <- ' #include <gsL/gsl_const_mksa.h> // decl of constants std::vector<double> volumes() { std::vector<double> v(5); v[0] = GSL_CONST_MKSA_US_GALLON; // 1 US gallon v[1] = GSL_CONST_MKSA_CANADIAN_GALLON; // 1 Canadian gallon v[2] = GSL_CONST_MKSA_UK_GALLON; // 1 UK gallon v[3] = GSL_CONST_MKSA_QUART; // 1 quart v[4] = GSL_CONST_MKSA_PINT; // 1 pint return v; }' gslVolumes <- cppFunction(code, depends="RcppGSL") mas aparece o seguinte erro file16e2b6cb966.cpp: In function ‘SEXPREC* sourceCpp_52966_volumes()’: file16e2b6cb966.cpp:30: error: ‘__result’ was not declared in this scope make: *** [file16e2b6cb966.o] Error 1 llvm-g++-4.2 -arch x86_64 -I/Library/Frameworks/R.framework/Resources/include -DNDEBUG -I/usr/local/include -I/Library/Frameworks/R.framework/Versions/3.0/Resources/library/Rcpp/include -I/usr/local/include -I"/Library/Frameworks/R.framework/Versions/3.0/Resources/library/Rcpp/include" -I"/Library/Frameworks/R.framework/Versions/3.0/Resources/library/RcppGSL/include" -fPIC -mtune=core2 -g -O2 -c file16e2b6cb966.cpp -o file16e2b6cb966.o Erro em sourceCpp(code = code, env = env, rebuild = rebuild, showOutput = showOutput, : Error 1 occurred building shared library.

1 1