Skip to content

cleanup aim1. minor file/ image changes. #122

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wrshoemaker wants to merge 2 commits into
base: master
Choose a base branch
from
Open

cleanup aim1. minor file/ image changes. #122

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
bd116af
cleanup aim1
wrshoemaker Mar 1, 2018
e8fed5f
Merge remote-tracking branch 'upstream/master'
wrshoemaker Mar 1, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
.Rhistory
.DS_Store
.Rapp.history
Aim1/.RData
Aim1/phylo/.RData
Binary file added Aim1/Rplot.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
17 changes: 10 additions & 7 deletions Aim1/code/test_weibull.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,8 @@ for(i in 1:length(strains)){
if (repObs["logabund"][[1]][2] - repObs["logabund"][[1]][1] > 1){
repObs <- repObs[-c(1), ]
}
repObs["prop"] <- repObs$logabund / repObs$logabund[1]
#repObs["prop"] <- repObs$logabund / repObs$logabund[1]
repObs["prop"] <- log(repObs$Abund / repObs$Abund[1])
#KBS0703.1$prop <- KBS0703.1$logabund/KBS0703.1$logabund[1]
# Initial parameters
#A = 200 # Initial death (larger = slower)
Expand All @@ -56,11 +57,11 @@ for(i in 1:length(strains)){
new.start[names(start) %in% names(mod.start)]<-mod.start
pscale<-as.numeric(new.start)
names(pscale)<-names(new.start)
#fit <- mle2(minuslogl=logabund ~ dnorm(mean = c * (time / a)^(b-1) * exp(-1*(time/a)^b), sd = z),
#fit <- mle2(minuslogl=prop ~ dnorm(mean = exp( -1 * ((time / a)^ b)), sd = z),
# start = new.start, data = repObs,
# control=list(parscale=pscale, maxit=1000),
# method="Nelder-Mead", hessian = T)
fit <- mle2(minuslogl=prop ~ dnorm(mean = exp( -1 * ((time / a)^ b)), sd = z),
# method="Nelder-Mead", hessian = T)
fit <- mle2(minuslogl=prop ~ dnorm(mean = -1 * ((time / a)^ b), sd = z),
start = new.start, data = repObs,
control=list(parscale=pscale, maxit=1000),
method="Nelder-Mead", hessian = T)
Expand Down Expand Up @@ -95,19 +96,21 @@ for(i in 1:length(strains)){

### *** Comment/Uncomment following code to make pdf figs*** ###
title=paste(strains[i]," rep ",reps[j])
plot(repObs$time,repObs$prop,main=title,ylim=c(0,1))
plot(repObs$time,repObs$prop,main=title,ylim=c(min(repObs$prop),0))
predTime=seq(0,max(repObs$time))
print(strains[i])
print(reps[j])
#exp( -1 * ((time / a)^ b))
#coef(best.fit)[3]
lines(repObs$time, exp( -1 * ((repObs$time / coef(best.fit)[1] )^ coef(best.fit)[2])),
#lines(repObs$time, exp( -1 * ((repObs$time / coef(best.fit)[1] )^ coef(best.fit)[2])),
# lwd=4, lty=2, col = "red")
lines(repObs$time, (-1 * ((repObs$time / coef(best.fit)[1] )^ coef(best.fit)[2])),
lwd=4, lty=2, col = "red")
counter=counter+1
}
}
}

dev.off()
summ=summ[!is.na(summ[,1]),]
#colnames(summ)=c('strain','rep','a','b','c','z','AIC', 'a.CI.2.5', 'a.CI.97.5', 'b.CI.2.5', 'b.CI.97.5', 'c.CI.2.5', 'c.CI.97.5', 'z.CI.2.5', 'z.CI.97.5')
Expand Down
123 changes: 123 additions & 0 deletions Aim1/data/weibull_log.R
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,123 @@
rm(list = ls())
getwd()
setwd("~/GitHub/LTDE/")

library('bbmle')

## Load Data
obs <- read.csv("data/longtermdormancy_20170620_nocomments.csv",
header = TRUE, stringsAsFactors = FALSE)
## Adding 1 to deal with log(0) observations
obs$Abund <- as.numeric(obs$Colonies) * 10 ^ as.numeric(obs$Dilution) + 1
strains <- sort(unique(obs$Strain))
strains <- strains[table(obs$Strain)>10]
obs <- obs[obs$Strain%in%strains,]
summ <- matrix(NA,length(strains)*max(obs$Rep),8)

pdf('output/decayFitsWeibull_log.pdf') # Uncomment to create pdf that will plot data and fits
counter <- 1

for(i in 1:length(strains)){
strainObs=obs[obs$Strain==strains[i],]

reps=unique(strainObs$Rep)
for(j in 1:length(reps)){

repObs=strainObs[strainObs$Rep==reps[j],]
# minimum of 10 data points
if(nrow(repObs)>10){
start=repObs[1,1]
time=(as.numeric(strptime(repObs$Firstread_date,format="%d-%b-%y",tz="EST"))-
as.numeric(strptime(start,format="%d-%b-%y",tz="EST")))/(3600*24)
repObs["time"] <- time + 1
repObs["logabund"] <- log10(repObs$Abund)
if (repObs["logabund"][[1]][2] - repObs["logabund"][[1]][1] > 1){
repObs <- repObs[-c(1), ]
}
#repObs["prop"] <- repObs$logabund / repObs$logabund[1]
repObs["prop"] <- log(repObs$Abund / repObs$Abund[1])
#KBS0703.1$prop <- KBS0703.1$logabund/KBS0703.1$logabund[1]
# Initial parameters
#A = 200 # Initial death (larger = slower)
#B = 1 # Bend (upper = 1 = first-order decay)
#C = round(max(repObs$logabund),1) # intercept
#Z = 6 # Error
grids<-list(a=c(1,10,50,100,200),b=c(0.1,0.5,1,1.1,1.5),z=c(0.1,1,10))
#start<-list(a=NA,b=NA,c=round(max(repObs$logabund),1),z=NA)
start<-list(a=NA,b=NA,z=NA)
grid.starts<-as.matrix(expand.grid(grids))
ncombos<-dim(grid.starts)[[1]]
# cycle through each combo
res.mat<-matrix(NA,nrow=ncombos,ncol=I(length(start)+1))
res.mod<-list()
for(k in 1:dim(grid.starts)[[1]]){
#some how need to match grid parameters to start lists.
mod.start<-as.list(grid.starts[k,])
new.start<-start
new.start[names(start) %in% names(mod.start)]<-mod.start
pscale<-as.numeric(new.start)
names(pscale)<-names(new.start)
#fit <- mle2(minuslogl=prop ~ dnorm(mean = exp( -1 * ((time / a)^ b)), sd = z),
# start = new.start, data = repObs,
# control=list(parscale=pscale, maxit=1000),
# method="Nelder-Mead", hessian = T)
fit <- mle2(minuslogl=prop ~ dnorm(mean = -1 * ((time / a)^ b), sd = z),
start = new.start, data = repObs,
control=list(parscale=pscale, maxit=1000),
method="Nelder-Mead", hessian = T)
res.mat[k,]<-c(coef(fit),AIC(fit))
res.mod[[k]]<-fit
}
colnames(res.mat)<-c(names(coef(fit)),"AIC")
best.fit<-res.mod[[which(res.mat[,'AIC']==min(res.mat[,'AIC']))[1]]]
summ[counter,1]=strains[i]
summ[counter,2]=reps[j]
#CIs <- confint( profile(best.fit))
print(coef(best.fit))
# a
summ[counter,3]=coef(best.fit)[1]
# b
summ[counter,4]=coef(best.fit)[2]
# c
summ[counter,5]=coef(best.fit)[3]
#summ[counter,5]= round(max(repObs$logabund),1)
# z
summ[counter,6]=AIC(best.fit)
#summ[counter,8]=CIs[1,1]
#summ[counter,9]=CIs[1,2]
#summ[counter,10]=CIs[2,1]
#summ[counter,11]=CIs[2,2]
#summ[counter,12]=CIs[3,1]
#summ[counter,13]=CIs[3,2]
#summ[counter,14]=CIs[4,1]
#summ[counter,15]=CIs[4,2]
summ[counter,7]=length(repObs$time)
ln_half <- log(0.5, base = exp(1)) *-1
half_life <- (ln_half ^ (1 / coef(best.fit)[2])) * coef(best.fit)[1]
summ[counter,8]=half_life


### *** Comment/Uncomment following code to make pdf figs*** ###
title=paste(strains[i]," rep ",reps[j])
plot(repObs$time,repObs$prop,main=title,ylim=c(min(repObs$prop),0))
predTime=seq(0,max(repObs$time))
print(strains[i])
print(reps[j])
#exp( -1 * ((time / a)^ b))
#coef(best.fit)[3]
#lines(repObs$time, exp( -1 * ((repObs$time / coef(best.fit)[1] )^ coef(best.fit)[2])),
# lwd=4, lty=2, col = "red")
lines(repObs$time, (-1 * ((repObs$time / coef(best.fit)[1] )^ coef(best.fit)[2])),
lwd=4, lty=2, col = "red")
counter=counter+1
}
}
}

dev.off()
summ=summ[!is.na(summ[,1]),]
#colnames(summ)=c('strain','rep','a','b','c','z','AIC', 'a.CI.2.5', 'a.CI.97.5', 'b.CI.2.5', 'b.CI.97.5', 'c.CI.2.5', 'c.CI.97.5', 'z.CI.2.5', 'z.CI.97.5')
colnames(summ)=c('strain','rep','a','b','z','AIC', 'N.obs', 'Half_life')

write.csv(summ,"data/weibull_log_results.csv")

101 changes: 101 additions & 0 deletions Aim1/data/weibull_log_results.csv
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,101 @@
"","strain","rep","a","b","z","AIC","N.obs"
"1","ATCC13985","1","57.5760113048932","0.67648262511964","0.897515049149791","179.029681892667","66"
"2","ATCC13985","2","50.8228410333533","0.644201373307315","0.908467873538774","180.631220734182","66"
"3","ATCC13985","3","3.17305345368234","0.358214889962649","1.28869850214531","193.330613847646","56"
"4","ATCC13985","4","8.05002131657765","0.416766253171727","0.936877799008328","157.612138836215","56"
"5","ATCC43928","1","5.08906729412192","0.423031576343636","1.29607036126235","227.536125949181","66"
"6","ATCC43928","2","8.58625316715676","0.471471439711458","1.43533808812812","241.00188147868","66"
"7","ATCC43928","3","1.41485560887317","0.274698753957306","0.906659488046098","153.942036224567","56"
"8","ATCC43928","4","124.531136267263","0.61162981709174","0.63938075276621","112.887169853393","55"
"9","KBS0701","1","0.0449413390487306","0.263012782350021","1.59460988662874","266.205099874135","69"
"10","KBS0701","2","0.159984974121919","0.287216902145944","1.68635022819535","273.931259764115","69"
"11","KBS0701","3","1.0497594698141","0.425249977737922","0.369343056791873","15.3045985044272","11"
"12","KBS0701","4","1.2442463430051","0.428792654871101","0.402569057401359","17.1990865310636","11"
"13","KBS0701","5","0.131806931731806","0.278002890499983","1.73975893289357","234.829508914809","58"
"14","KBS0701","6","0.11861522072844","0.288229408528618","1.38242796183582","211.652228659638","59"
"15","KBS0702","1","10.553061858749","0.531021706239314","0.886779451930214","138.471095581888","51"
"16","KBS0702","2","8.63697988475462","0.484014297796612","0.818017013717281","130.235123611459","51"
"17","KBS0702","3","7.18568453924163","0.412938197742558","0.773732446759374","119.913191841697","49"
"18","KBS0702","5","11.7517966580184","0.511033137343314","0.827330270820816","126.48584877759","49"
"19","KBS0702","4","8.67805425772814","0.433604295014221","0.889622223624492","136.197603591598","50"
"20","KBS0702","6","12.5026190162885","0.542028937060911","0.867514675541646","131.127998422537","49"
"21","KBS0703","1","3.09590436947271","0.340640962042937","0.714292175919755","116.413938253049","51"
"22","KBS0703","2","204.766867290655","0.983815886588752","0.90260509298486","140.278397610749","51"
"23","KBS0703","3","5.86542863003184","0.367030540005129","0.579480566359455","91.5877599999787","49"
"24","KBS0703","4","2.64128753116559","0.324309456067429","0.554024147869495","85.5265562184569","48"
"25","KBS0703","6","9.57779857995625","0.494935083542353","0.728062784509772","67.6841899839353","28"
"26","KBS0703","5","16.0838616321163","0.609374351916796","1.01653680862307","80.6397039324249","26"
"27","KBS0705","1","4.74293340675482","0.462629407642157","1.45384312108225","110.003261504252","29"
"28","KBS0705","2","0.14301245319875","0.208737808871396","1.47482707035669","179.515452476575","48"
"29","KBS0705","3","0.994434362080044","0.285733178401649","1.28939176413509","166.623445691713","48"
"30","KBS0705","4","0.777819705693486","0.31085472636745","1.69573897136951","192.913763250491","48"
"31","KBS0706","1","64.401877329416","0.59196061276618","0.972336666669765","133.966777135542","46"
"32","KBS0706","2","111.566319964714","0.75558899307325","1.04691124318224","140.760496902622","46"
"33","KBS0706","3","78.8381985003132","0.590582253689269","0.81550377956972","117.776489755705","46"
"34","KBS0706","4","113.557092333617","0.682821757759137","0.781061830990009","113.810595069442","46"
"35","KBS0707","1","2.20449411126757","0.276444252010417","0.908301021487707","132.985227615991","48"
"36","KBS0707","2","0.864035529616693","0.241059774825411","1.3594109136175","171.689956640154","48"
"37","KBS0707","3","0.913921265519429","0.264327465494098","1.33184199263796","169.727506261229","48"
"38","KBS0707","4","2.60556407852199","0.279396602689294","1.15183114458258","155.789004893529","48"
"39","KBS0710","1","1.39539707213863","0.348517214847981","1.13608788936692","157.564509395764","49"
"40","KBS0710","2","0.565255295937808","0.304271114806702","1.1184040053062","156.019920953613","49"
"41","KBS0710","3","3.31862882607408","0.428474247364105","1.11038064903934","158.364621874005","50"
"42","KBS0710","4","2.79107311619321","0.436279090045552","1.17964805780078","164.412224768936","50"
"43","KBS0711","1","0.0333446528581724","0.239753118511919","1.09636647599496","226.595772850905","73"
"44","KBS0711","2","0.252372226571199","0.273616425539473","0.871392248667837","193.067626389763","73"
"45","KBS0711","3","0.0269210605257621","0.221667646394692","0.898840533283133","197.600146386777","73"
"46","KBS0711","4","1.46694691251382","0.441858574019229","0.551157852271984","29.0507158776858","14"
"47","KBS0711","11","303.198612763851","0.807699916846199","0.643109581569513","125.260793350742","61"
"48","KBS0711","10","546.28953899869","2.03834109997829","0.799504164255993","135.083998562092","54"
"49","KBS0711W","1","0.157721774938032","0.2585075244385","0.997804596827395","164.670394416954","56"
"50","KBS0711W","2","0.643183763467108","0.328696398493759","1.40388011215741","202.917590253268","56"
"51","KBS0711W","3","0.440418730735792","0.282198609844244","1.18797977626173","184.211379203736","56"
"52","KBS0711W","4","0.394288448598924","0.263995087365449","1.21874530206582","187.072799422251","56"
"53","KBS0712","1","1.26093445944646","0.316012870709948","1.13735925850481","151.478416394877","47"
"54","KBS0712","2","0.122857478046521","0.218395991563584","0.891873807873325","131.224887519365","48"
"55","KBS0712","3","0.337441538317373","0.251460764381184","0.916658229298621","131.200278146358","47"
"56","KBS0712","4","0.00910828560500057","0.157219972698626","1.01587958289539","143.734150364831","48"
"57","KBS0713","3","0.014815723271025","0.179563372025191","1.16731507008911","163.359893464607","50"
"58","KBS0713","4","0.0146679071226549","0.194864350942151","1.01983824978665","149.853076511322","50"
"59","KBS0713","1","0.0171763130544269","0.193365183651568","0.940847563889355","139.079823065734","49"
"60","KBS0713","2","0.109236606107578","0.240021682159323","0.849469821243195","129.068473260819","49"
"61","KBS0714","1","15.7031456236137","0.895927080762207","1.55169702438518","95.1956695750156","24"
"62","KBS0714","2","17.3497350985668","0.920857986883279","1.49478807207806","93.4039714748551","24"
"63","KBS0714","3","13.5581063827781","0.842355527950507","1.50524012696145","93.7378683893821","24"
"64","KBS0714","4","15.892136938892","0.905806886525179","1.56500015715344","95.6068868668269","24"
"65","KBS0715","1","0.000168733457962035","0.144662174884253","1.03617534135487","209.625538510555","70"
"66","KBS0715","2","0.000122111955379719","0.144324860383312","1.02225419129159","207.735725879127","70"
"67","KBS0715","3","5.74360853675006e-05","0.137914114340581","1.08495626780349","216.067128971644","70"
"68","KBS0715","4","0.000321860782762532","0.153906876709854","1.01887613653882","207.269921406325","70"
"69","KBS0721","1","0.812031532614833","0.318572308860505","0.982516789457468","193.774841183716","67"
"70","KBS0721","2","0.0841743441752453","0.24240160597916","1.2605961927689","227.164582981485","67"
"71","KBS0721","3","0.0115126004155363","0.209373394501207","1.26029272902989","227.135669419342","67"
"72","KBS0721","4","0.171264891451397","0.265117240079671","1.34960767553099","236.314149363414","67"
"73","KBS0722","1","6.97200907944259","0.373203905951593","0.940677170517635","136.343782333159","48"
"74","KBS0722","2","1.77028099520983","0.266951114559207","0.966653219991871","138.964508131064","48"
"75","KBS0722","3","2.93058727906981","0.302330158878131","1.0367157503133","145.681515872988","48"
"76","KBS0722","4","2.20415222366699","0.301725244735672","1.19214679395231","155.901918323045","47"
"77","KBS0724","1","8.95703517463506","0.489957354733044","1.04024046437729","143.083404412554","47"
"78","KBS0724","2","8.17783656610518","0.473974737399186","0.897517816228768","129.219799498386","47"
"79","KBS0724","3","12.4538294752803","0.541552768383641","0.888018526649462","128.216648903807","47"
"80","KBS0724","4","7.05530468487364","0.473931378184008","1.02896822601535","142.063746916409","47"
"81","KBS0725","1","12.1418678942695","0.674024199105972","1.3484189901913","164.042214747553","46"
"82","KBS0725","2","18.4931963167587","0.761454161045611","1.30700410342094","167.92601582475","48"
"83","KBS0725","3","14.0594773951975","0.723898988018977","1.44310157520879","173.854103369557","47"
"84","KBS0725","4","12.787801663052","0.709803798875661","1.52293157583265","182.599285166462","48"
"85","KBS0727","1","1.09000979964027","0.370965029560497","1.1593231091073","115.673672542194","35"
"86","KBS0727","2","6.42929131706983","0.549059774995095","1.534406739466","164.853040590854","43"
"87","KBS0727","3","15.2803632440957","0.700132202227439","1.55308617833511","206.794868860032","54"
"88","KBS0727","4","11.7408272585699","0.651917900831077","1.47548602451037","157.863991030497","42"
"89","KBS0801","1","29.7871542864078","0.570561309988683","0.737167457124751","104.031480738548","44"
"90","KBS0801","2","25.2763669986262","0.556695269061455","0.714669096416174","101.305617272967","44"
"91","KBS0801","3","25.7926393860597","0.552683274659608","0.913459414453118","122.905659629987","44"
"92","KBS0801","4","63.8108738947654","0.693004534867913","1.0278008755366","130.391115222109","43"
"93","KBS0802","1","5.39610015149686","0.421948262366496","1.62790665527564","257.630085480631","66"
"94","KBS0802","2","4.53349150417438","0.403696629259194","1.75262067269996","267.363918495918","66"
"95","KBS0802","3","6.67404044411606","0.43959634009492","1.64865822790301","259.299993057265","66"
"96","KBS0802","4","6.0377456602531","0.43293604745193","1.77044853901695","264.718465470923","65"
"97","KBS0812","1","371926.632251891","0.143800735488994","0.326394479501541","46.1062213207067","67"
"98","KBS0812","2","26551.2156282938","0.603353439035324","0.324391243171935","45.286166557052","67"
"99","KBS0812","3","24238.5156067062","0.52957739209423","0.373588000643459","64.2057887219245","67"
"100","KBS0812","4","9263991475.28429","0.100759662238973","0.412880192360696","77.6234726622928","67"
Loading