###Endereço : ## ' http://leg.ufpr.br/~paulojus/embrapa/Rembrapa/Rembrapa26.html' ## Carregando as funções que iram integrar o pacote source('gamfit.R') source('hch2d.R') source('waddmod.R') source('kerreg2d.R') source('cvkreg2d.R') source('gamfit.sig.R') source('gam.sim.R') source('gam.sig.start.R') source('inout.R') source('npts.R') source('areapl.R') source('my.plotrisk.R') source('my.polymap.R') source('plotrisk.sil.R') source('polymap.sil.R') source('saveplot.R') source('image.legend.R') ### Criando o .tar.gz package.skeleton(name='spatial',list=c('cvkreg2d','gamfit','hch2d','kerreg2d','waddmod','gam.sig.start','gam.sim','gamfit.sig','inout','npts','areapl','saveplot','polymap.sil','plotrisk.sil','my.polymap','my.plotrisk','image.legend')) # construindo o pacote ### No terminal do linux R CMD build 'nome do pacote' ### Agora só no R ## Instalando o pacote no R install.packages('spatial_1.0.tar.gz',repos=NULL) ###Carregando os bancos source('dump.casos.neo.new') source('dump.casos.pos') source('dump.contpoa') source('dumpcontrol') ###Carregando o pacote require(spatial) require(splancs) ##Ajustando o GAM #1. Crie um objetgamfit.neo<-gamfit(form=form,neo,h=0,hvals=hvals,ngrid=4000)o contendo as coordenadas dos casos e controles pts<-rbind(cbind(casos.neo.new$cx,casos.neo.new$cy),cbind(control$cx,control$cy)) #2. Crie um objeto contendo o indicador de casos e controles y <- c(rep(1,nrow(casos.neo.new)),rep(0,nrow(control))) #3. Escolha os valores de banda candidatos para a validacao cruzada hvals <- seq(100, 500, length = 15) #4. Criando matrix de dados #Covariaveis para os casos ob.var<-data.frame(sexo=casos.neo.new$SEXO,peso=casos.neo.new$PESO,idademae=casos.neo.new$IDADMAE,instrmae=casos.neo.new$INSTRMAE,ges=as.numeric(casos.neo.new$ges),tipograv=casos.neo.new$TIPOGRAV,tipopart=as.numeric(casos.neo.new$tipopart1)) co.var<-data.frame(sexo=control$SEXO,peso=control$PESO,idademae=control$IDADMAE,instrmae=control$INSTRMAE,ges=as.numeric(control$ges),tipograv=control$TIPOGRAV,tipopart=as.numeric(control$tipopart1)) #Unindo covariaveis de casos e controles num data.frame var<-rbind(ob.var,co.var) gam.neo<-data.frame(y,var) # Formula para ajuste do modelo no formato usual do R form<-(y~sexo+peso+idademae+instrmae+ges+tipograv+tipopart) gamfit.neo<-gamfit(form=form,gam.neo,pts=pts,region=contpoa,h=0,hvals=hvals,ngrid=4000) ls() #Para visualizar os efeitos fixos estimados gamfit.neo$beta #Para visualizar o efeito espacial estimado digite polymap(contpoa) image(gamfit.neo$g2est,add=TRUE) polymap(contpoa,add=TRUE) #Repita os procedimentos acima para ajustar GAM as casos pos ########################################################################### # Obtendo os contornos de tolerancia e o teste global da hipotese nula de # risco nao variavel espacialmente ########################################################################### start<-gam.sig.start(form,gam.neo,pts,region=contpoa,h=gamfit.neo$h,ngrid=4000) #primeiras 10 simulacoes sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=NULL) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) #mais 10 simulacoes sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) sim<-gam.sim(form,gam.neo,pts,region=contpoa,start,ngrid=4000,prev.sims=sim) #acrescente mais simulacoes se desejar #idealmente 99 para um nivel de 5% e 499 para um nivel de 1%) neo.sig<-gamfit.sig(gamfit.neo,region=contpoa,sim) my.plotrisk(neo.sig,legend=TRUE) #pvalor do teste global de risco contante na regiao neo.sig$pval # se o teste for significativo entao podemos interpretar os contornos calculados #######TESTANDO NOS DADOS DE CURITIBA########################################## ############################################################################### neo <- read.table('neo_contr1.txt',h=T) head(neo) names(neo) pts <- -rbind(cbind(neo$XCOORD[neo$Y==1],neo$YCOORD[neo$Y==1]),cbind(neo$XCOORD[neo$Y==0],neo$YCOORD[neo$Y==0])) y <- factor(neo$Y) hvals <- seq(100, 500, length = 25) #4. Criando matrix de dados #Covariaveis para os casos ob.var<-data.frame(sexo=neo$SEXO[neo$Y==1],peso=neo$PESO[neo$Y==1],idademae=neo$IDADEMAE[neo$Y==1],instrmae=neo$ESCMAE[neo$Y==1],ges=neo$GESTACAO[neo$Y==1],tipograv=neo$GRAVIDEZ[neo$Y==1],tipopart=neo$PARTO[neo$Y==1]) #Covariaveis para os controles co.var<-data.frame(sexo=neo$SEXO[neo$Y==0],peso=neo$PESO[neo$Y==0],idademae=neo$IDADEMAE[neo$Y==0],instrmae=neo$ESCMAE[neo$Y==0],ges=neo$GESTACAO[neo$Y==0],tipograv=neo$GRAVIDEZ[neo$Y==0],tipopart=neo$PARTO[neo$Y==0]) #Unindo covariáveis de casos e controles num data.frame var<-rbind(ob.var,co.var) gam.neo<-data.frame(y,var) # Formula para ajuste do modelo no formato usual do R form<-"y~sexo+peso+idademae+instrmae+ges+tipograv+tipopart" #Ajustando GAM com todas as covariaveis gamfit.neo<-gamfit(form=form,gam.neo,pts=pts,h=0,hvals=hvals,ngrid=4000)