Commit ff637ffe authored by Marta Pittavino's avatar Marta Pittavino
Browse files

=DiB materials: data, codes and figures

parent dc53a66b
##################################################################################################
# R file which peforms 2560 bootstrap analysis for cohort sheep, lepto pomona data!!!
##################################################################################################
setwd("~/Desktop/ABN_GLM_Comparison_DataAndCodes/R/LeptoPomonaJags") ### Set the proper working directory
library( abn);
library( coda);
load("~/Desktop/ABN_GLM_Comparison_DataAndCodes/RData/leptopom_noethn.RData")
orig.data <- mydatpomdef_noethn;
max.par <- 7;
start <- seq( 1,2560,by=80)
stop <- seq( 80,2560,by=80);
mydists<-list( Lepto ="binomial",
Sex = "binomial",
Hunt = "binomial",
Farm = "binomial",
Kill = "binomial",
Glov = "binomial",
Glass = "binomial",
Mask = "binomial",
Age = "gaussian",
Time = "gaussian",
Work1 = "binomial",
Work2 = "binomial",
Work3 = "binomial",
Plant1 = "binomial",
Plant2 = "binomial",
Plant3 = "binomial",
Plant4= "binomial"
);
ban <- matrix( 0,17,17);
colnames( ban) <- rownames( ban) <- names( orig.data);
ban["Work1","Work2"]<-1;
ban["Work1","Work3"]<-1;
ban["Work2","Work1"]<-1;
ban["Work2","Work3"]<-1;
ban["Work3","Work1"]<-1;
ban["Work3","Work2"]<-1;
ban["Plant1","Plant2"]<-1;
ban["Plant1","Plant3"]<-1;
ban["Plant2","Plant1"]<-1;
ban["Plant2","Plant3"]<-1;
ban["Plant3","Plant1"]<-1; # now ban arc from work1 to work3
ban["Plant3","Plant2"]<-1; # now ban arc from work2 to work3
ban["Plant1","Plant4"]<-1; # now ban arc from to
ban["Plant4","Plant1"]<-1; # now ban arc from to
ban["Plant2","Plant4"]<-1;
ban["Plant4","Plant2"]<-1;
ban["Plant3","Plant4"]<-1;
ban["Plant4","Plant3"]<-1;
ban["Plant3","Work1"]<-1;
ban["Work1","Plant3"]<-1;
retain <- matrix( 0,17,17);
colnames( retain)<-rownames( retain)<-names( orig.data);
dags<-list();
#########################################
for(i in start[index]:stop[index]){
########################################
#Create bootstrap data
#1. create parameter file with unique random seed
init.file<-paste("init_",i,sep="");
cat(paste("\".RNG.name\" <-\"base::Mersenne-Twister\"","\n",sep=""),file=init.file,append=FALSE);
cat(paste("\".RNG.seed\" <- ",i,"\n",sep=""),file=init.file,append=TRUE);
#2. create script file with unique output file name
run.file<-paste("script_",i,sep="");
cat("model in LeptoPomonaModel_NoEthn.bug
data in leptopomona_postparams_noethn.R
compile, nchains(1)
",file=run.file);
cat(paste("parameters in ",init.file,"\n",sep=""),file=run.file,append=TRUE);
cat("initialize
update 10000
monitor Lepto, thin(10)
monitor Sex, thin(10)
monitor Hunt, thin(10)
monitor Farm, thin(10)
monitor Kill, thin(10)
monitor Glov, thin(10)
monitor Glass, thin(10)
monitor Mask, thin(10)
monitor Age, thin(10)
monitor Time, thin(10)
monitor Work1, thin(10)
monitor Work2, thin(10)
monitor Work3, thin(10)
monitor Plant1, thin(10)
monitor Plant2, thin(10)
monitor Plant3, thin(10)
monitor Plant4, thin(10)
update 3840, by(1000)
",file=run.file,append=TRUE);
out.file<-paste("out_",i,sep="");
cat(paste("coda *, stem(\"",out.file,"\")\n",sep=""),file=run.file,append=TRUE);
#3. run the MCMC sampler
system(paste("jags ",run.file,sep=""));
#4. read in mcmc data and convert to format suitable for mostprobable
boot.data<-read.coda(paste(out.file,"chain1.txt",sep=""),paste(out.file,"index.txt",sep=""));
boot.data<-as.data.frame(boot.data);
for(j in 1:dim(orig.data)[2]){if(is.factor(orig.data[,j])){boot.data[,j]<-as.factor(boot.data[,j]);
levels(boot.data[,j])<-levels(orig.data[,j]);}}
#5. run the MostProb search on the bootstrap data
boot1.cache<-buildscorecache( data.df=boot.data,data.dists=mydists, max.parents=max.par,
dag.banned=ban, dag.retained=retain, verbose=F);
dags[[i]]<-mostprobable( score.cache=boot1.cache);
unlink(c(init.file,run.file,out.file,paste(out.file,"chain1.txt",sep=""),paste(out.file,"index.txt",sep="")));#tidy up
datadir <- '~/Desktop/ABN_GLM_Comparison_DataAndCodes/RData/noethn/Bootstrapping/'
save( dags, file=paste( datadir,"mp2560boot_noethn",index,".RData",sep=""));
}
\ No newline at end of file
#!/bin/bash
# 2560Bootstrapping_NoEthn.R
Rout="Rout/noethn"
Rfiles="R/LeptoPomonaJags"
echo
# Standard syntax.
for a in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
do
rm -f $Rfiles/"torun.r"
echo "index <- $a" >> $Rfiles/torun.r
cat $Rfiles/2560Bootstrapping_NoEthn.R >> $Rfiles/torun.r
R CMD BATCH $Rfiles/torun.r $Rout/out2560boot$a.txt
echo -n "$a"
done
rm -f $Rfiles/"torun.r"
echo; echo
rm(list=ls())
setwd("~/Desktop/ABN_GLM_Comparison_DataAndCodes") ### Set the proper working directory
load("RData/leptopom_noethn.RData")
library(abn);
mydata <- mydatpomdef_noethn;
mydists <- list( Lepto ="binomial",
Sex = "binomial",
Hunt = "binomial",
Farm = "binomial",
Kill = "binomial",
Glov = "binomial",
Glass = "binomial",
Mask = "binomial",
Age = "gaussian",
Time = "gaussian",
Work1 = "binomial",
Work2 = "binomial",
Work3 = "binomial",
Plant1 = "binomial",
Plant2 = "binomial",
Plant3 = "binomial",
Plant4= "binomial"
);
### fit model GLOBAL BEST SCORE DAG - this is mp.dag_noethn7: mp.dag_noethn
# Lepto Sex Hunt Farm Kill Glov Glass Mask Age Time Work1 Work2 Work3 Plant1 Plant2 Plant3 Plant4
# Lepto 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0
# Sex 0 0 1 0 0 0 0 0 0 1 1 0 1 0 1 1 1
# Hunt 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0
# Farm 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# Kill 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0
# Glov 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# Glass 0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 1 0
# Mask 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0
# Age 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
# Time 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0
# Work1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# Work2 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0
# Work3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# Plant1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0
# Plant2 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
# Plant3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# Plant4 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
mydag <- matrix( c( ### Matrix of DAG Definition. Rows are children, Cols are parents.
# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
#Lepto Sex Hunt Farm Kill Glov Glass Mask Age Time Work1 Work2 Work3 Plant1 Plant2 Plant3 Plant4
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, # Lepto
0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, # Sex
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, # Hunt
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # Farm
0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, # Kill
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # Glov
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, # Glass
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, # Mask
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, # Age
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, # Time
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # Work1
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, # Work2
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # Work3
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # Plant1
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, # Plant2
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # Plant3
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 # Plant4
),byrow=TRUE,ncol=17);
colnames( mydag)<-rownames( mydag)<-names( mydata);
## use fitabn just to check mydag is correct (no typos as mlik should = -3655.078
print( fitabn(dag.m=mydag,data.df=mydata,data.dists=mydists)$mlik);## ok. have correct model
#################################################################################################
### now compute marginals for each and every parameter in the model
#################################################################################################
marg.f<-fitabn( dag.m=mydag,data.df=mydata,data.dists=mydists,
compute.fixed=TRUE,n.grid=1000);
library(Cairo);
CairoPDF("MargPlots_LeptoPomona_NoEthn.pdf");
for(i in 1:length(marg.f$marginals)){
cat("processing marginals for node:",nom1<-names(marg.f$marginals)[i],"\n");
cur.node<-marg.f$marginals[i];
cur.node<-cur.node[[1]];
for(j in 1:length(cur.node)){
cat("processing parameter:",nom2<-names(cur.node)[j],"\n");
cur.param<-cur.node[[j]];
plot(cur.param,type="l",main=paste(nom1,":",nom2));
}
}
dev.off();
marnew <- marg.f$marginals[[1]]
for(i in 2: length(marg.f$marginals)){
marnew <- c(marnew, marg.f$marginals[[i]])
}
myarea<-rep(NA,length(marnew));names(myarea)<-names(marnew);
for(i in 1:length(marnew)){
tmp<-spline(marnew[[i]]);
myarea[i]<-sum(diff(tmp$x)*tmp$y[-1]);
}
library(Cairo);
mycols<-rep("green",length(marnew));
mycols[1:4]<-"red";# Lepto
mycols[5:12]<-"orange"; # Sex
mycols[13:14]<-"yellow"; # Hunt
mycols[15]<-"blue"; # Farm
mycols[16:19] <-"lightblue"; # Kill
mycols[20]<-"mistyrose"; # Glov
mycols[21:26]<-"lightcyan"; # Glass
mycols[27:30]<- "lavender"; # Mask
mycols[31:33]<-"cornsilk"; # Age
mycols[34:37]<-"brown" ; # Time
mycols[38]<-"darkgreen"; # Work1
mycols[39:41] <-"cyan"; # Work2
mycols[42]<-"purple"; # Work3
mycols[43:44]<-"darkblue"; # Plant1
mycols[45:46]<- "darkorange"; # Plant2
mycols[47]<-"darkred"; # Plant3
CairoPNG("Area_LeptoPomona_NoEthn.png",pointsize=10,width=720,height=640);
par(las=2);
par(mar=c(8.1,4.1,4.1,2.1));
barplot(myarea,ylab="Area under Density",ylim=c(0,1.5), col=mycols);
dev.off();
margs <- marnew;
mymat<-matrix(rep(NA,length(margs)*3),ncol=3);rownames(mymat)<-names(margs);colnames(mymat)<-c("2.5%","50%","97.5%");
ignore.me<-union(grep("\\(Int",names(margs)),grep("prec",names(margs)));
comment<-rep("",length(margs));
for(i in 1:length(margs)){
tmp<-margs[[i]];
tmp2<-cumsum(tmp[,2])/sum(tmp[,2]);
mymat[i,]<-c(tmp[which(tmp2>0.025)[1]-1,1],
tmp[which(tmp2>0.5)[1],1],
tmp[which(tmp2>0.975)[1],1]);
myvec<-mymat[i,];
if( !(i%in%ignore.me) && (myvec[1]<0 && myvec[3]>0)){comment[i]<-"not sig. at 5%";}
mymat[i,]<-as.numeric(formatC(mymat[i,],digits=3,format="f"));
}
mymat
### Now create the DATA for going to JAGS for the MCMC: Important point!!!
print(names( marnew));
m <- marnew;
Lepto.p <- cbind( m[["Lepto|(Intercept)"]], m[["Lepto|Work1"]], m[["Lepto|Work2"]], m[["Lepto|Work3"]]);# Lepto --> 3 parents;
Sex.p <- cbind( m[["Sex|(Intercept)"]], m[["Sex|Hunt"]], m[["Sex|Time"]], m[["Sex|Work1"]],
m[["Sex|Work3"]], m[["Sex|Plant2"]], m[["Sex|Plant3"]], m[["Sex|Plant4"]]); # Sex --> 7 parents;
Hunt.p <- cbind( m[["Hunt|(Intercept)"]], m[["Hunt|Age"]]); #Hunt --> 1 parent;
Farm.p <- cbind( m[["Farm|(Intercept)"]]); # Farm --> NO PARENT, Orfan Node! Changed wrt Initial Model.
Kill.p <- cbind( m[["Kill|(Intercept)"]], m[["Kill|Sex"]], m[["Kill|Farm"]], m[["Kill|Mask"]]); # Kill --> 3 parents;
Glov.p <- cbind( m[["Glov|(Intercept)"]]); # Glov --> NO PARENT, Orfan Node!
Glass.p <- cbind( m[["Glass|(Intercept)"]], m[["Glass|Work1"]], m[["Glass|Work2"]],
m[["Glass|Work3"]], m[["Glass|Plant2"]], m[["Glass|Plant3"]]); # Glass --> 5 parents;
Mask.p <- cbind( m[["Mask|(Intercept)"]], m[["Mask|Work1"]],
m[["Mask|Work2"]], m[["Mask|Work3"]]); # Mask --> 3 parents;
Age.p <- cbind(m[["Age|(Intercept)"]], m[["Age|Time"]]); # Age --> 1 parents: Gaussian;
prec.Age.p <- m[["Age|precision"]];
Time.p <- cbind(m[["Time|(Intercept)"]], m[["Time|Work3"]], m[["Time|Plant2"]]); # Time --> 2 parents: Gaussian;
prec.Time.p <- m[["Time|precision"]];
Work1.p <- cbind( m[["Work1|(Intercept)"]]); # Work1 --> NO PARENT, Orfan Node!
Work2.p <- cbind( m[["Work2|(Intercept)"]], m[["Work2|Glov"]], m[["Work2|Plant3"]]); # Work2 --> 2 parents!
Work3.p <- cbind( m[["Work3|(Intercept)"]]); # Work3 --> NO PARENT, Orfan Node!
Plant1.p <- cbind( m[["Plant1|(Intercept)"]], m[["Plant1|Glass"]]); # Plant1 --> 1 parent!
Plant2.p <- cbind( m[["Plant2|(Intercept)"]], m[["Plant2|Mask"]]); # Plant2 --> 1 parent!
Plant3.p <- cbind( m[["Plant3|(Intercept)"]]); # Plant3 --> NO PARENT, Orfan Node!
Plant4.p <- cbind( m[["Plant4|(Intercept)"]], m[["Plant4|Mask"]]); # Plant4 --> 1 parent!
datadir <- 'R/LeptoPomonaJags/'
dump(c("Lepto.p", "Sex.p", "Hunt.p", "Farm.p", "Kill.p", "Glov.p", "Glass.p", "Mask.p",
"Age.p", "prec.Age.p", "Time.p", "prec.Time.p", "Work1.p", "Work2.p",
"Work3.p", "Plant1.p", "Plant2.p", "Plant3.p", "Plant4.p"),
file=paste( datadir, "leptopomona_postparams_noethn.R", sep=""));
proc.time()
#user system elapsed
#60.749 34.262 67737.084
rm(list=ls())
setwd("~/Desktop/ABN_GLM_Comparison_DataAndCodes") ### Set the proper working directory
load("RData/leptopom_noethn.RData")
mydata <- mydatpomdef_noethn;
library(abn);
N <- 2560;
############################################################################################
mydag <- matrix( c( ### Matrix of DAG Definition. Rows are children, Cols are parents.
# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
#Lepto Sex Hunt Farm Kill Glov Glass Mask Age Time Work1 Work2 Work3 Plant1 Plant2 Plant3 Plant4
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, # Lepto
0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, # Sex
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, # Hunt
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # Farm
0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, # Kill
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # Glov
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, # Glass
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, # Mask
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, # Age
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, # Time
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # Work1
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, # Work2
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # Work3
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # Plant1
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, # Plant2
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, # Plant3
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 # Plant4
),byrow=TRUE,ncol=17);
colnames( mydag)<-rownames( mydag)<-names( mydata);
mydists <- list( Lepto ="binomial",
Sex = "binomial",
Hunt = "binomial",
Farm = "binomial",
Kill = "binomial",
Glov = "binomial",
Glass = "binomial",
Mask = "binomial",
Age = "gaussian",
Time = "gaussian",
Work1 = "binomial",
Work2 = "binomial",
Work3 = "binomial",
Plant1 = "binomial",
Plant2 = "binomial",
Plant3 = "binomial",
Plant4= "binomial"
);
## use fitabn just to check mydag is correct (no typos as mlik should = -3655.078)
print( fitabn(dag.m=mydag,data.df=mydata,data.dists=mydists)$mlik);## ok. have correct model
sum( mydag); # 30 arcs, check notes
bestdag <- mydag;
#############################################################################################
## read ALL files with mp[number].RData and create one list of results.:
boot.dags<-list();
these<-grep("mp2560boot_noethn\\d+.RData", dir());
num<-1;
for(i in dir()[these]){
load(i);
tmp<-dags[which(unlist(lapply(dags,sum))>0)];
for(j in 1:length(tmp)){
boot.dags[[num]]<-tmp[[j]];num<-num+1;}
rm(dags);
}
if(FALSE){
scores<-rep(0,length(boot.dags));for(i in 1:length(boot.dags)){scores[i]<-fitabn(dag.m=boot.dags[[i]],data.df=mydata,data.dists=mydists)$mlik;}
scores.b<-scores[-which(scores< -N)];
orig.score<-fitabn(dag.m=bestdag,data.df=mydata,data.dists=mydists)$mlik;
plot(density(scores.b,from=min(scores.b),to=max(scores.b)))
abline(v=orig.score,lwd=2,col="blue")
}
boot.dags.trim<-boot.dags;
for(i in 1:length(boot.dags)){
boot.dags.trim[[i]]<-boot.dags.trim[[i]]*bestdag;
}
arc.freq<-lapply(boot.dags.trim,sum);arc.freq<-table(unlist(arc.freq));
library(Cairo);
CairoPNG("LeptoFreqBootRes2K.png",pointsize=10,width=720,height=700);
par(las=1);
par(mar=c(6.1,6.1,4.1,2.1));
barplot(arc.freq,ylab="",xlab="",col="skyblue",names.arg=names(arc.freq),ylim=c(0,400),cex.axis=1.5);
par(las=1);
#axis(1,at=seq(1,22,by=1));
mtext("Number of arcs in bootstrap DAG",1,line=3,cex=1.5);
par(las=3);
mtext("Frequency out of 2000",2,line=4,cex=1.5)
dev.off();
arc.freq<-lapply(boot.dags.trim,sum);arc.freq<-table(unlist(arc.freq));
par(mai=c(1.1,1.4,0.15,0.05))
par(las=1);
barplot(arc.freq,ylab="",xlab="",col="black",names.arg=names(arc.freq),ylim=c(0,350),cex.axis=1.1);
par(las=1);
mtext("Number of arcs in bootstrap DAG",1,line=3,cex=1.5);
par(las=3);
mtext("Frequency out of 2560",2,line=4,cex=1.5)
dev.off();
##############################################################################################################################
total.dag<-matrix(rep(0,dim(bestdag)[2]^2),ncol=dim(bestdag)[2]);colnames(total.dag)<-rownames(total.dag)<-colnames(bestdag);
for(i in 1:length(boot.dags)){
if(sum(boot.dags[[i]])>0){total.dag<-total.dag+boot.dags[[i]];}}
total.dag<-total.dag*bestdag;
total.dag
f<-function(val,limit){if(val<limit){return(0);} else {return(1);}}
bestdag.trim<-apply(total.dag,c(1,2),FUN=f,limit=N/2);
bestdag.trim.nodir<-bestdag;
bestdag.trim.nodir[,]<-0;
child<-NULL;parent<-NULL;
for(i in 1:dim(total.dag)[1]){
for(j in 1:dim(total.dag)[2]){
if(i>j){
if(total.dag[i,j]>total.dag[j,i]){m.i<-i;m.j<-j;} else {m.i<-j;m.j<-i;}
if(total.dag[i,j]+total.dag[j,i]>N/2){
bestdag.trim.nodir[m.i,m.j]<-1;}
}}}
all.equal(bestdag.trim, bestdag.trim.nodir)
tographviz(dag.m=bestdag.trim,data.df=mydata,data.dists=mydists,outfile="BootLeptoPomonaDAG2K.dot");
system("dot -Tpng -o BootLeptoPomonaDAG2K.png BootLeptoPomonaDAG2K.dot");
system("evince BootLeptoPomonaDAG2K.png&")
save( bestdag.trim,file=paste("Bestdag_leptopomona_trim2K.RData",sep=''));
\ No newline at end of file
rm(list=ls())
setwd("~/Desktop/ABN_GLM_Comparison_DataAndCodes") ### Set the proper working directory
load("RData/leptopom_noethn.RData")
load("RData/noethn/Bootstrapping/Bestdag_leptopomona_trim2K.RData");
mydata <- mydatpomdef_noethn;
library(abn);
mydists <- list( Lepto ="binomial",
Sex = "binomial",
Hunt = "binomial",
Farm = "binomial",
Kill = "binomial",
Glov = "binomial",
Glass = "binomial",
Mask = "binomial",
Age = "gaussian",
Time = "gaussian",
Work1 = "binomial",
Work2 = "binomial",
Work3 = "binomial",
Plant1 = "binomial",
Plant2 = "binomial",
Plant3 = "binomial",
Plant4= "binomial"
);
marg.f<-fitabn(dag.m=bestdag.trim,data.df=mydata,data.dists=mydists,compute.fixed=TRUE,n.grid=1000);
library(Cairo);
CairoPDF("LeptoPomona_PostBootPlots2K.pdf");
for(i in 1:length(marg.f$marginals)){
cat("processing marginals for node:",nom1<-names(marg.f$marginals)[i],"\n");
cur.node<-marg.f$marginals[i];## get marginal for current node - this is a matrix [x,f(x)]
cur.node<-cur.node[[1]];# this is always [[1]] for models without random effects
for(j in 1:length(cur.node)){
cat("processing parameter:",nom2<-names(cur.node)[j],"\n");
cur.param<-cur.node[[j]];
plot(cur.param,type="l",main=paste(nom1,":",nom2));
}
}
dev.off();
marnew <- marg.f$marginals[[1]]
for(i in 2: length(marg.f$marginals)){
marnew <- c(marnew, marg.f$marginals[[i]])
}
margs <- marnew;
mymat <- matrix(rep(NA,length(margs)*3),ncol=3);rownames(mymat)<-names(margs);colnames(mymat)<-c("2.5%","50%","97.5%");
ignore.me <- union(grep("\\(Int",names(margs)),grep("prec",names(margs)))
comment <- rep("",length(margs));
for(i in 1:length(margs)){
tmp <- margs[[i]];
tmp2 <- cumsum(tmp[,2])/sum(tmp[,2]);
mymat[i,] <- c(tmp[which(tmp2>0.025)[1]-1,1],
tmp[which(tmp2>0.5)[1],1],
tmp[which(tmp2>0.975)[1],1]);
myvec <- mymat[i,];
if( !(i%in%ignore.me) && (myvec[1]<0 && myvec[3]>0)){comment[i]<-"not sig. at 5%";}
mymat[i,] <- as.numeric(formatC(mymat[i,],digits=3,format="f"));
}
cbind( mymat);