oddsAndEnds/GDAnalysis.R at master · ewong347/oddsAndEnds · GitHub

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
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
#Analyze Growth data from tn93GD.R

### USAGE: Rscript GDAnalysis.R tn93GDOutput.RData ###

##Stat Presentation functions
#__________________________________________________________________________________________#

#Obtain the GAIC, a measure of fit between predicted and actual growth
gaic <- function(res=res)  {
  stat <- sapply(1:ncol(res), function(x) {
    #Extract full and fit data
    fit <- res[[1,x]]
    full <- res[[2,x]]

    #Place growth and forecast data in dfs for fit and full growth
    #df1 <- data.frame(Growth = fit$growth, Pred = fit$forecast)
    #df2 <- data.frame(Growth = full$growth, Pred = full$forecast)

    #Place growth and forecast data in dfs for fit and full growth
    df1 <- data.frame(Growth = fit$growth, Pred = fit$forecast)
    df2 <- data.frame(Growth = fit$growth, Pred = fit$csize * (sum(fit$growth)/sum(fit$csize)))

    #Model growth as a function of forecast for fit and full growth models
    mod1 <- glm(Growth ~ Pred, data = df1, family = "poisson")
    mod2 <- glm(Growth ~ Pred, data = df2, family = "poisson")

    #Calculate GAIC
    mod1$aic-mod2$aic
  })

  #Present and return data
  plot(colnames(res), stat, ylab = "GAIC", xlab = "Cutoff", type= "l", cex.lab=1.5)
  return(stat)
}

#__________________________________________________________________________________________#

#Expecting the output from a run of tn93GD.Rdata
args = commandArgs(trailingOnly = T)

#Loads the output from tn93GD.RData
load(args)

## Skeleton Code For plotting
stat <- sapply(1:ncol(res), function(x) {
  fit <- res[[1,x]]

})
plot(colnames(res), stat, ylab = "" , xlab= "")