text:
######################################
### process year by year estimates
######################################

processEsts <- function(ot,ageDist18){
    for(i in 1:length(ot)){
    ### add year as data
        for(j in 1:length(ot[[i]])){
            ot[[i]][[j]]$year <- i
        }
    ### separate race/sex estimates by sex
    for(lev in c('hs','cc','bach'))
        for(sex in 1:2){
            ot[[i]][[paste0(lev,'Race',c('M','F')[sex])]] <-
                subset(ot[[i]][[paste0(lev,'RaceSex')]],SEX==sex,select=-SEX)
        }
    }

    ## remname se column
    for(i in 1:length(ot)){
        for(j in 1:length(ot[[i]])){
            names(ot[[i]][[j]]) <- gsub('se2','se',names(ot[[i]][[j]]))
            if('raceEth'%in%names(ot[[i]][[j]])){
                ot[[i]][[j]]$raceEth <- gsub('Hispanic','Latinx',ot[[i]][[j]]$raceEth)
                #ot[[i]][[j]]$raceEth <- gsub(' ',' ',ot[[i]][[j]]$raceEth)
            }
        }
    }

    if(!missing(ageDist18)) attr(ot,'weightDat') <- ageDist18

    ot
}

#### combine separate year-by-year sets of estimates into
### estimates from one set of subgroups
combineDat <- function(nnn,ot,deaf=NULL){
    tdat <- do.call('rbind', lapply(ot, function(x) x[[nnn]]))
    names(tdat) <- gsub('bach|hs','y',names(tdat))
    names(tdat)[grep('edLev',names(tdat))] <- 'y'

    for(i in 1:ncol(tdat)) if(is.factor(tdat[,i])) tdat[,i] <- as.character(tdat[,i])

    if('SEX'%in%names(tdat)) tdat$SEX <- c('Male','Female')[tdat$SEX]

    if(!is.null(deaf)){
        tdat <- tdat[tdat$DEAR==ifelse(deaf,1,2),]
        tdat$DEAR <- NULL
    }

    tdat$se[tdat$se==0] <- NA

    subCols <- setdiff(names(tdat),c('y','AGEP','se','Freq','se1','se2','year','ww'))
    subsets <- if(length(subCols)>1) do.call(expand.grid,lapply(tdat[,subCols],unique))
               else cbind(unique(tdat[,subCols]))

   if(!is.null(attr(ot,'weightDat'))){
        tdat <- full_join(tdat,attr(ot,'weightDat'))
    } else
        tdat$ww <- if(any(is.na(tdat$se))) tdat$Freq else 1/tdat$se^2

    list(tdat=tdat,subCols=subCols,subsets=subsets)
}

##########################################
### functions to make tables
##########################################

###################
## functions for attainment rates
###################
getRates <- function(nnn,ot,years=c('08','18'),digits=2){
    rnd <- function(x) sprintf(paste0('%.',digits,'f'),x)
    tab <- map_dfr(years,function(y)
        ot[[y]][[nnn]]%>%
            mutate(
                rate=paste0(rnd(`I(attain >= edLev)TRUE`),' (',rnd(se),')'),
                year=paste0('20',y)
                )%>%
            select(-`I(attain >= edLev)TRUE`,-se,-Freq)
                   )%>%
        tidyr::pivot_wider(names_from=year,values_from=rate)
    if('SEX'%in%names(tab)) tab$SEX <- c('Male','Female')[tab$SEX]
    if(endsWith(nnn,'Tot')){
        tab <- mutate(tab,subgroup='Overall')%>%select(DEAR,subgroup,everything())
    } else names(tab)[2] <- 'subgroup'

    tab
}

ratesTab <- function(lev,subs,ot,years,digits,onlyDeaf){
    tab <- map_dfr(subs,function(ss)
                       getRates(paste0(lev,ss),ot=ot,years=years,digits=digits)%>%
                       filter(!subgroup%in%c('35-44','45-54','55-64'))%>%
                       mutate(SUB=ss)
                   )
    if(onlyDeaf) tab <- filter(tab,DEAR==1)%>%select(-DEAR)
    tab
}

###################
## functions for attainment trends
###################

#### makes trend table given combineDat output
trends <- function(tdat, subCols){
    trnds <- tdat%>%
        group_by(!!! lapply(subCols,sym))%>%
            group_modify(estTrend1)

    trnds$subgroup <- if(length(subCols)==1) "Overall" else trnds[[subCols[2]]]

    trnds
}

### estimate the trend
estTrend1 <- function(ddd,...){
    form <- y~year+as.factor(AGEP)
    trnd <- summary(lm_robust(form,data=ddd,weights=ww))$coef['year',]
    if(is.nan(trnd['Std. Error']))
        trnd <- summary(lm_robust(form,data=ddd,weights=ww,se_type='HC1'))$coef['year',]
    as.data.frame(rbind(trnd))
}

trendRegressionTab <- function(ddd,...){
  form <- y~year+as.factor(AGEP)
  trnd <- lm_robust(form,data=ddd,weights=ww)
trnd
}

###################
## functions for gap trends
###################

### estimate one gap trend
gapTrendOne <- function(tdat,...){
  mod <-
    if(any(is.na(tdat$se))){
      lm_robust(y~(DEAR+splines::ns(AGEP,5))*year,tdat,weights=ww)
    } else
      lm_robust(y~(DEAR+as.factor(AGEP))*year,weights=ww,data=tdat)

  as.data.frame(summary(mod)$coef)['DEAR:year',c(1,2,4)]
}

### estimate gap trend for one set of subgroups, given combineDat output
gapTrend <- function(tdat,subCols,...){
    gaps <-
        if(length(subCols)==1){
            gapTrendOne(tdat)
        } else
            tdat%>%
                group_by_at(vars(!! sym(subCols[2])))%>%
                    group_modify(gapTrendOne)
    gaps$subgroup <- if(length(subCols)==1) "Overall" else gaps[[subCols[2]]]
    gaps
}

###################
## functions to compute all tables and put them together
###################

### regression output -> estimate (se)
tabClean <- function(.data,name,digits=2){
    rnd <- function(x) sprintf(paste0('%.',digits,'f'),x)
    if(!'DEAR'%in%names(.data)) .data$DEAR <- 1
    .data%>%
       ungroup()%>%
       transmute(
           DEAR,
           subgroup,
           !! name:=paste0(rnd(Estimate),' (',rnd(`Std. Error`),')'),
           !! paste0(name,'p') := `Pr(>|t|)`
           )#%>%dplyr::select(SUB,everything())
}

### makes table for trend or gap for one subgroup/level combination
trendTab1 <- function(lev,SUB,ot,trendOrGap=c('trend','gap'),digits=2){
    fun <- switch(trendOrGap,trend=trends,gap=gapTrend)
    out <- combineDat(paste0(lev,SUB),ot)%$%
        fun(tdat,subCols)%>%
            tabClean(trendOrGap,digits=digits)
    #if(SUB=='RaceF') out$subgroup <- paste(out$subgroup,'F')
    #if(SUB=='RaceM') out$subgroup <- paste(out$subgroup,'M')
    if(SUB=='ByAgeCat') out <- filter(out,subgroup=='25-34')
    out$SUB <- SUB
    out
}

### p-value stars
stars <- function(pval)
    ifelse(pval<0.001,'***',
           ifelse(pval<0.01,'**',
                  ifelse(pval<0.05,'*',
                         ifelse(pval<0.1,'.',''))))

### put together tables for all sets of subgroups, levels, attainment or gap.
### does multiplicity adjustment, returns list of tables for 3 levels
tabsFun <- function(ot, trendOrGap=c('trend','gap'),
                     subs=c('Tot','ByAgeCat','Sex','Race','RaceM','RaceF'),
                     onlyDeaf=TRUE,digits=2){

    if(trendOrGap=='gap') onlyDeaf <- TRUE

    tabs <- sapply(c('hs','cc','bach'),
                   function(lev){
                       trnds <- map_dfr(subs,~trendTab1(lev=lev,.,ot=ot,trendOrGap=trendOrGap))

                       if(onlyDeaf){
                           trnds <- trnds%>%filter(DEAR==1)%>%select(-DEAR)
                       } else
                           trnds <- trnds%>%
                               mutate(deaf=ifelse(DEAR==1,'Deaf','Hearing'))%>%
                                   select(-DEAR)%>%
                                       pivot_wider(names_from=deaf,values_from=c('trend','trendp'))%>%
                                           rename(trend=trend_Deaf,trendp=trendp_Deaf)
                       trnds
                   }, simplify=FALSE
                   )

    psTot <- map_dbl(tabs,~ .[[paste0(trendOrGap,'p')]][1])%>%p.adjust(.,'holm')
    psSub <- do.call('c',map(tabs,~ .[[paste0(trendOrGap,'p')]][-1]))%>%p.adjust(.,'fdr')

    for(lev in c('hs','cc','bach'))
        tabs[[lev]][[paste0(trendOrGap,'p')]] <- c(psTot[lev],psSub[startsWith(names(psSub),lev)])

    tabs <- map(tabs,function(x) x%>%mutate(!! sym(trendOrGap):=paste0(!!sym(trendOrGap),stars(!!sym(paste0(trendOrGap,'p')))))%>%
                    select(-starts_with('trendp'),-starts_with('gapp')))

    tabs

}

#### makes all the trend tables
allTabs <- function(overTimeAge,overTime,
                    subs=c('Tot','ByAgeCat','Sex','Race','RaceM','RaceF'),
                    onlyDeaf=TRUE,digits=2,years=c('08','18')){

    tabs <- sapply(c('trend','gap'),
                function(x) tabsFun(overTimeAge,trendOrGap=x,subs=subs,onlyDeaf=onlyDeaf,digits=digits),
                   simplify=FALSE)

    rates <- sapply(c('hs','cc','bach'),ratesTab,
                    subs=subs,ot=overTime,years=years,digits=digits,onlyDeaf=onlyDeaf,
                    simplify=FALSE)

    sapply(c('hs','cc','bach'),
           function(x)
              full_join(rates[[x]],tabs$trend[[x]])%>%
              full_join(tabs$gap[[x]])%>%
              rename(`Deaf\nAttainment`=trend,`Deaf-Hearing\nGap`=gap)%>%
              select(SUB,subgroup,everything())%>%
              arrange(match(SUB,subs))%>%
              mutate(SUB=c(Tot="",ByAgeCat="",Sex="Sex",Race="Race/\nEthnicity",RaceM="Race/\nEthnicity:\nMale",RaceF="Race/\nEthnicity:\nFemale")[SUB]),
           simplify=FALSE
           )
  }

##########################################
### functions to make plots
##########################################

### for each level, four plots: overall, 25-34, by sex, by race/eth

figFun <- function(lev,ot,subs=c('Tot','ByAgeCat','Sex','Race')){

    dd <- map_dfr(subs, function(ss){
                ccc <- combineDat(paste0(lev,ss),ot)
                tdat <- ccc$tdat
                if(ss=='Race') tdat$raceEth <-
                    c(`African American`="Afr. Am.", `American Indian`="Am. Ind.",  `Asian/PacIsl`="Asn/PacIsl", `Latinx`="Latinx",  `Other`="Other", `White`="White")[tdat$raceEth]
                if(ss%in%c('Tot','ByAgeCat')) tdat$subgroup=' ' else names(tdat)[2] <- 'subgroup'
                if(ss=='ByAgeCat') tdat <- subset(tdat,ageRange=='25-34')
                tdat$SUB <- switch(ss,
                                   Tot='Overall',
                                   ByAgeCat='Age 25-34',
                                   Sex='By Gender',
                                   Race='By Race/Ethnicity'
                                   )
                #tdat$lab <- ifelse(tdat$DEAR==1,tdat$subgroup,paste0(tdat$subgroup,' '))
                tdat$lab <- ifelse(tdat$year==11,tdat$subgroup,'')
                tdat$grp <- paste(tdat$DEAR,tdat$subgroup)
                tdat$deaf <- ifelse(tdat$DEAR==1,'Deaf','Hearing')
                tdat
            })#%>%
        dd$SUB=factor(dd$SUB,levels=unique(dd$SUB))#)%>%
        ggplot(dd,aes(year,y,linetype=deaf,group=grp,label=lab))+
            geom_line()+geom_point(size=0.5)+
                scale_x_continuous(
                    breaks=1:11,
                    minor_breaks=NULL,
                    labels=c(2008:2018),
                    limits=c(1,13)
                    )+
                #geom_dl(method=list("last.qp",cex=0.5))+
                geom_text_repel(size=2,force=0.25,min.segment.length=0.25,segment.size=0.5)+
                scale_y_continuous(labels=function(x) paste0(x,'%'))+
                facet_wrap(~SUB,ncol=2)+
                labs(
                  title=paste("Figure", switch(lev,hs=1,cc=2,bach=3)),
                  subtitle=(paste(switch(lev, hs='High School',cc='Associates',bach='Bachelors'),'Degree Attainment, 2008-2018')),
                    x=NULL,
                    y=paste('% Attaining', switch(lev, hs='High School',cc='Associates',bach='Bachelors'),'Degree or Higher'),
                    linetype=NULL
                    )+
                    theme(legend.position='bottom',legend.margin=margin(t=0,b=-.3,unit='cm'),
                        axis.text.x = element_text(angle = 90,hjust = 1, vjust = 0.5),
                          plot.title = element_text(face = "bold"),
                          plot.subtitle=element_text(face="italic"))
}

##########################################
### results text functions
##########################################

############
### results table
############
pfun <- function(p)
  ifelse(p>=0.01,paste0("p=",sprintf("%.2f",p)),
         ifelse(p>0.001,"p<0.01","p<0.001"))

tab1 <- function(lev,SUB,what,ot=if(what=='lev') overTime else overTimeAge,deaf=what!="gap",digits=2,tab=FALSE){

  fun <- switch(what,trend=trends,gap=gapTrend,lev=levs)
  out <- combineDat(paste0(lev,SUB),ot)%$%
    fun(tdat,subCols)
  if(deaf) out <- filter(out,DEAR==1)

  if(tab) return(out)

  if("Estimate"%in%names(out))
    out <- out%>%
      arrange(desc(Estimate))%>%
        mutate(Estimate=round(Estimate,2))
  if("y2008" %in% names(out))
    out <- mutate_at(out,vars(starts_with("y")),~paste0(round(.,digits=1),'%'))

  if(SUB=='ByAgeCat') out <- filter(out,ageRange=='25-34')

  if("Pr(>|t|)"%in%names(out)) out <- mutate(out,pval=pfun(`Pr(>|t|)`))

  out <- as.data.frame(out)

  if(n_distinct(out[,2])==nrow(out)) rownames(out) <- out[,2]

  out
}

levs <- function(tdat,subCols){
  tdat%>%
    filter(year%in%c(1,11))%>%
    mutate(year=2007+year)%>%
    select(!!!syms(subCols),y,se,year)%>%
    pivot_wider(names_from='year',values_from=c("y","se"),names_sep='')
}

interactionTest <- function(out){
  w <- 1/out$`Std. Error`^2
  T <- sum(w*out$Estimate)/sum(w)
  Q <- sum(w*(out$Estimate-T)^2)
  pchisq(Q,length(w)-1,lower.tail=FALSE)
}

pval <- function(top,bottom)
    2*pnorm(as.numeric(-abs(top/bottom)))

diffP <- function(ests,ses){
    est <- ests[2]-ests[1]
    se <- sqrt(ses[2]^2+ses[1]^2)
    pval(est,se)
  }

pairWise <- function(out,groupName=names(out)[2]){
  out <- out%>%arrange(desc(Estimate))
  pvals <- NULL
  for(i in 1:(nrow(out)-1))
    for(j in (i+1):nrow(out))
      pvals <- c(pvals,
                 setNames(
                   diffP(out$Estimate[c(i,j)],out$`Std. Error`[c(i,j)]),
                   paste0(out[i,groupName],'>',out[j,groupName])
                   )
                 )
  #p.adjust(pvals,'hommel')
  pvals
}

totTest1 <- function(lev,SUB,what,ot=overTimeAge){
  if(SUB%in%c('Sex','Race')){
    out <- tab1(lev=lev,SUB=SUB,what=what,ot=ot,tab=TRUE)
    if(nrow(out)==1) return(NULL)
    if(nrow(out)==2) pvals <- c(sex=diffP(out$Estimate,out$`Std. Error`)) else
    pvals <- c(overall=interactionTest(out),pairWise(out))
  } else
    pvals <- bind_rows(
      tab1(lev,'RaceM',what)%>%mutate(sex=paste0(raceEth,'M')),
      tab1(lev,'RaceF',what)%>%mutate(sex=paste0(raceEth,'F'))
    )%>%
    group_by(raceEth)%>%group_map(~pairWise(.,groupName='sex'))%>%do.call('c',.)
  data.frame(levs=lev,subs=SUB,what=what,test=names(pvals),pvals=pvals)
}

totTest <- function(overTimeAge,levs=c('hs','cc','bach'),subs=c('Sex','Race','RaceSex'),what=c('trend','gap')){
  tab <- expand.grid(levs=levs,subs=subs,what=what,stringsAsFactors=FALSE)
  pvals <- map_dfr(1:nrow(tab),~totTest1(tab$levs[.],tab$subs[.],tab$what[.],overTimeAge))
  pvals$pvals[pvals$test=='overall'] <- p.adjust(pvals$pvals[pvals$test=='overall'],'holm')

  pvals <- pvals%>%
    group_by(levs,subs,what)%>%
    mutate(
      overall=ifelse('overall'%in%test,pvals[test=='overall']<0.05,TRUE)
      )%>%
      ungroup()
  pvals$overall[pvals$test=='overall'] <- FALSE
  pvals$pvals[pvals$overall] <- p.adjust(pvals$pvals[pvals$overall],'fdr')
  pvals$pvals[!pvals$overall&pvals$test!='overall'] <- 1
  pvals
}

extrSubs <- function(ccc,varb){
    ## tab <- with(ccc,{
    ##          tab <- NULL
    ##          for(i in 1:nrow(subsets)){
    ##              keep <- rep(TRUE,nrow(tdat))
    ##              for(j in 1:ncol(subsets)){
    ##                  keep <- keep&(tdat[,subCols[j]]==subsets[i,j])
    ##              }
    ##              tab <- cbind(tab,tdat[[varb]][keep][order(tdat$year[keep])])
    ##          }
    ##          tab})
    ddd <- ccc$tdat[,c('year',varb,ccc$subCols)]
    names(ddd) <- c('year','varb',paste0('sub',seq(length(ccc$subCols))))
    tab <- spread(ddd, key=sub1,value=varb)
    rownames(tab) <- 2007+tab$year
    tab$year <- NULL
    tab
}

tabFunBegin <- function(nnn,ot,deaf=NULL){
    ccc <- combineDat(nnn,ot,deaf)

    trnd <- trends(nnn,ot)
    if(!is.null(deaf)) trnd <- trnd[grep(paste0('DEAR=',ifelse(deaf,1,2)),names(trnd))]

    if(nnn%in%names(overTimeAge[[1]])){
        trndsAdj <- trends(nnn,overTimeAge)
        if(!is.null(deaf)) trndsAdj <-
            trndsAdj[grep(paste0('DEAR=',ifelse(deaf,1,2)),names(trndsAdj))]
    }

    ests <- extrSubs(ccc,'y')*1
    ses <- extrSubs(ccc,'se')*1

    ests <- rbind(ests,
                  Growth=ests[nrow(ests),]-ests[1,])

    ses <- rbind(ses,Growth=sqrt(ses[nrow(ses),]^2+ses[1,]^2))

    out <- list(ests=ests,ses=ses,trnd=trnd[sapply(names(ests),function(nn) grep(nn,names(trnd)))])
    if(nnn%in%names(overTimeAge[[1]]))
        out$trndsAdj <- trndsAdj[sapply(names(ests),function(nn) grep(nn,names(trndsAdj)))]
    out
}

estsRS <- function(tabb){
    tab <- rbind(tabb$ests,Trend=NA,`Trend (Adj)`=NA)
    ses <- rbind(tabb$ses,Trend=NA,`Trend (Adj)`=NA)
    for(j in 1:ncol(tab)){
        nm <- colnames(tab)[j]
        sex <- ifelse(grepl('Male',nm),'Male','Female')
        race <- substr(nm,nchar(sex)+2,nchar(nm))
        trd <- tabb$trnd[[grep(paste0('SEX=',sex,' raceEth=',race,' DEAR='),names(tabb$trnd))]]
        trdA <- tabb$trndsAdj[[grep(paste0('SEX=',sex,' raceEth=',race,' DEAR='),names(tabb$trndsAdj))]]
        tab['Trend',j] <- trd[1]*1
        tab['Trend (Adj)',j] <- trdA[1]*1
        ses['Trend',j] <- trd[2]*1
        ses['Trend (Adj)',j] <- trdA[2]*1
    }
    list(ests=tab,ses=ses)
}

tabFun <- function(nnn,ot,deaf=NULL,moe=TRUE,retEst=TRUE){

    tfb <- tabFunBegin(nnn,ot,deaf)
    if(grepl('RaceSex',nnn)) tfb <- estsRS(tfb)
    ests <- tfb$ests
    ses <- tfb$ses
    if(!grepl('RaceSex',nnn)){
        trnd <- tfb$trnd
        trndsAdj <- tfb$trndsAdj

        ests <- rbind(ests,
                      Trend=vapply(1:ncol(ests),function(j)
                          trnd[[grep(colnames(ests)[j],names(trnd))]][1]*1,1),
                      `Trend (Adj)`=vapply(1:ncol(ests),function(j)
                          trndsAdj[[grep(colnames(ests)[j],names(trnd))]][1]*1,1))
        ses <- rbind(ses,
                     Trend=vapply(1:ncol(ests),function(j)
                         trnd[[grep(colnames(ests)[j],names(trnd))]][2]*1,1),
                     `Trend (Adj)`=vapply(1:ncol(ests),function(j)
                         trndsAdj[[grep(colnames(ests)[j],names(trnd))]][2]*1,1))
    }
  #    if(retEst)
  return(list(ests=ests,ses=ses))

  #  tabFinish(ests,ses,moe)
}

tabFinish <- function(anal1,deaf=NULL,moe=TRUE){ #ests,ses,moe,strs){
    if(is.null(deaf)){
        ests <- anal1$ests
        ses <- anal1$ses
        strs <- cbind(anal1$trendStars,c('',''))
        rownames(strs) <- c('unadj','adj')
        diffStars <- anal1$diffStars
        ss <- if('SS'%in%names(anal1)) anal1$SS else c(anal1$deafSS,anal1$hearSS)
    } else if(deaf){
        ests <- anal1$deaf$ests
        ses <- anal1$deaf$ses
        strs <- anal1$trendStars
        ss <- anal1$deafSS
    } else{
        ests <- anal1$hear$ests
        ses <- anal1$hear$ses
        strs <- rbind(unadj=rep('',ncol(ests)),adj=rep('',ncol(ests)))
        ss <- anal1$hearSS
    }
    tab <- ests
    for(j in 1:ncol(ests)){
        for(i in 1:nrow(ests)){
            if(moe) tab[i,j] <- paste(round(ests[i,j],1),'$\\pm$',round(ses[i,j]*1.96,1))
            else tab[i,j] <- paste0(round(ests[i,j],1),' (',round(ses[i,j],1),')')
        }
        tab['Trend',j] <- paste0(tab['Trend',j],strs['unadj',j])
        tab['Trend (Adj)',j] <- paste0(tab['Trend (Adj)',j],strs['adj',j])

    }
    tab <- rbind(tab,`Sample Size/year`= format(ss,big.mark=','))
    if(is.null(deaf)){
        colnames(tab) <- c('Deaf','Hearing')
        ## tab <- cbind(tab,Difference=
        ##                    c(
        ##                      vapply(1:nrow(ests),
        ##                        function(i)
        ##                          diffFun(ests[i,],ses[i,],moe,
        ##                            ifelse(rownames(ests)[i]=='Trend (Adj)',diffStars[2],'')),
        ##                        'a'
        ##                      ),''
        ##                    )
        ## )
    }

  ## c(rep('',(nrow(tab)-4)),
  ##                            diffFun(ests['Growth',],ses['Growth',],moe,''),
  ##                            diffFun(ests['Trend',],ses['Trend',],moe,diffStars[1]),
  ##                            diffFun(ests['Trend (Adj)',],ses['Trend (Adj)',],moe,diffStars[2]),''))
  ##   }

    tab <- tab%>%
      tibble::rownames_to_column('Year')

    tab
}

diffFun <- function(ests,ses,moe,diffStars=''){
    est <- ests[1]-ests[2]
    se <- sqrt(ses[2]^2+ses[1]^2)
    p <- pval(est,se)

    out <- round(est,1)
    out <- if(moe) paste(out,'$\\pm$',round(1.96*se,1)) else paste0(out,' (',round(se,1),')')
    out <- paste0(out,diffStars)
    out
}

    ## errbar(2008:2016,estsD[,1],estsD[,1]+2*sesD[,1],estsD[,1]-2*sesD[,1],
##            ylim=c(0,1),type='l')#,...)
##     lines(2008:2016,estsH[,1],lty=3)
##     errbar(2008:2016,estsH[,1],estsH[,1]+2*sesH[,1],estsH[,1]-2*sesH[,1],lty=2,add=TRUE)
##     for(i in 1:ncol(estsD)){
##         errbar(2008:2016,estsD[,i],estsD[,i]+2*sesD[,i],estsD[,i]-2*sesD[,i],col=i,add=TRUE,type='l')
##         lines(2008:2016,estsH[,i],col=i,lty=3)
##         errbar(2008:2016,estsH[,i],estsH[,i]+2*sesH[,i],estsH[,i]-2*sesH[,i],col=i,add=TRUE,lty=2)
##     }
##     legend(pos,legend=c(colnames(estsD),'','Deaf','Hearing'),col=c(1:ncol(estsD),0,1,1),lty=c(rep(1,ncol(estsD)),0,1,2))
## }

changeFig <- function(nnn,ot,...){
    cccD <- combineDat(nnn,ot,deaf=TRUE)
    cccH <- combineDat(nnn,ot,deaf=FALSE)

    estsD <- extrSubs(cccD,'y')
    sesD <- extrSubs(cccD,'se')
    estsH <- extrSubs(cccH,'y')
    sesH <- extrSubs(cccH,'se')

    for(i in nrow(estsD):1){
        estsD[i,] <- 100*(estsD[i,]-estsD[1,])
        estsH[i,] <- 100*(estsH[i,]-estsH[1,])
    }

    ests <- rbind(melt(estsD),melt(estsH))
    names(ests) <- c('Year','sub','Change')
    ests$deaf <- c(rep('Deaf',prod(dim(estsD))),rep('Hearing',prod(dim(estsH))))

    ggplot(ests,aes(Year,Change,color=sub,shape=deaf,linetype=deaf))+
        geom_point()+geom_line()+labs(color='',shape='',linetype='')+
        geom_hline(yintercept=0,linetype=2)+
        ylab('Change Since 2008 (Percentage Points)')+xlab(NULL)+
        theme(axis.text.x = element_text(angle = 90, hjust = 1,vjust=0.5))

}

##     lll <- range(c(unlist(estsD,estsH)))
##     plot(2008:2016,estsD[,1],pch=16,ylim=lll,type='b',...)
##     points(2008:2016,estsH[,1],pch=17,lty=2,type='b')
##     for(i in 1:ncol(estsD)){
##         points(2008:2016,estsD[,i],pch=16,col=i,type='b')
##         points(2008:2016,estsH[,i],col=i,pch=17,lty=2,type='b')
##         #deafBeta <- trendsUnadj[[nnn]][[paste0(cccD$subCols,'=',colnames(estsD)[i],' DEAR=1')]][1]
##         #hearBeta <- trendsUnadj[[nnn]][[paste0(cccD$subCols,'=',colnames(estsD)[i],' DEAR=2')]][1]

##         #abline(a=-2008*deafBeta,b=deafBeta,col=i)
##         #abline(a=-2008*hearBeta,b=hearBeta,col=i,lty=2)
##         abline(h=0,lty=2)
##     }
##     legend(pos,legend=c(colnames(estsD),'','Deaf','Hearing'),col=c(1:ncol(estsD),0,1,1),lty=c(rep(1,ncol(estsD)),0,1,2),pch=c(rep(NA,ncol(estsD)+1),16,17))
## }

resultsTab <- function(ttt,fname='trends.csv'){
    out <- data.frame(lev=character(0),
                      deaf=character(0),
                      sss=character(0),
                      est=numeric(0),
                      serr=numeric(0),
                      pp=character(0),
                      stars=character(0))
    for(nnn in names(ttt)){
        rrr <- ttt[[nnn]]
        lev <- rep(ifelse(grepl('hs',nnn),'High School','Bachelors'),length(rrr))
        deaf <- ifelse(grepl('DEAR=1',names(rrr)),'Deaf','Hearing')
        sss <- gsub('DEAR=[12]','',names(rrr))
        if(grepl('25',nnn)) sss <- paste0(sss,'ageRange=25-29')
        rrr <- do.call('rbind',rrr)
        ddd <- data.frame(lev,deaf,sss,round(rrr[,1:2]*1,round(abs(log(min(rrr[,2]),10)))-2),
                          ifelse(rrr[,4]<0.001,'<0.001',round(rrr[,4],3)),
                          ifelse(rrr[,4]<0.001,'***',ifelse(rrr[,4]<0.01,'**',
                                                            ifelse(rrr[,4]<0.05,'*',
                                                                   ifelse(rrr[,4]<0.1,'.','')))))

        names(ddd) <- names(out)
        out <- rbind(out,ddd)
    }
    write.table(out,paste0('output/',fname),sep=',',row.names=FALSE,
                col.names=c('Level','Deaf or Hearing','Sub-Group',
                    'Change per Year','Std. Error','p-value',''))
    invisible(out)
}

#estsD,estsH,sesD,sesH,moe=TRUE,diffStars){
diffDeafHear <- function(anal1,moe=TRUE){
    estsD <- anal1$deaf$ests
    estsH <- anal1$hear$ests
    sesD <- anal1$deaf$ses
    sesH <- anal1$hear$ses

    diffTab <- NULL
    diffStars <- rbind('',anal1$diffStars)
    rownames(diffStars) <- c('Growth','Trend','Trend (Adj)')
    for(varb in c('Growth','Trend','Trend (Adj)')){
        diff <- NULL
        for(j in 1:ncol(estsD)){
            diff <- c(diff,diffFun(c(estsD[varb,j],estsH[varb,j]),
                                   c(sesD[varb,j],sesH[varb,j]),moe,diffStars[varb,j]))
        }
        diffTab <- rbind(diffTab,diff)
    }
    rownames(diffTab) <- c('Growth','Trend','Trend (Adj)')
    colnames(diffTab) <- colnames(estsD)
    diffTab
}

diffSub <- function(ests,ses,varb,moe=TRUE){
    diffTab <- matrix('',ncol(ests)-1,ncol(ests)-1)

    for(i in 1:(ncol(ests)-1)){
        for(j in (i+1):ncol(ests)){
            diffTab[i,j-1] <- diffFun(c(ests[varb,i],ests[varb,j]),
                                    c(ses[varb,i],ses[varb,j]),moe)
        }
    }
    rownames(diffTab) <- colnames(ests)[-ncol(ests)]
    colnames(diffTab) <- colnames(ests)[-1]
    diffTab
}

toList <- function(lst,vec){
  fac <- factor(
    do.call(
      'c',
      lapply(
        names(lst),
        function(nm) if(is.matrix(lst[[nm]])) rep(nm,nrow(lst[[nm]])) else nm
      )
    )
  )
    split(vec,fac)
}

mult1 <- function(trend,diff,alpha,adj){
    colnum <- ifelse(adj,2,1)
    trendTot <- trend[[1]][colnum]
    diffTot <- diff[[1]][colnum]
    trend2 <- lapply(trend[-1],function(x) x[,colnum])
    diff2 <- lapply(diff[-1],function(x) x[,colnum])

    if(trendTot>alpha) return(list(trend=lapply(trend,function(x) rep(FALSE,length(x))),
                                   diff=lapply(trend,function(x) rep(FALSE,length(x)))))

    newPs <- p.adjust(do.call('c',trend2),method='hommel')
    rej <- c(TRUE,newPs<=alpha)

    diffPs <- c(diffTot,do.call('c',diff2))
    rejD <- rep(FALSE,length(diffPs))
    rejD[rej] <- p.adjust(diffPs[rej],method='fdr')<=alpha

    list(trend=toList(trend,rej),diff=toList(diff,rejD))
}

mult2 <- function(pTrendL,pDiffL,alpha,adj){
    pHS <- pTrendL[grep('hs',names(pTrendL))]
    diffHS <- pDiffL[grep('hs',names(pTrendL))]

    pCC <- pTrendL[grep('cc',names(pTrendL))]
    diffCC <- pDiffL[grep('cc',names(pTrendL))]

    pBach <- pTrendL[grep('bach',names(pTrendL),ignore.case=TRUE)]
    diffBach <- pDiffL[grep('bach',names(pTrendL),ignore.case=TRUE)]

    rejHS <- mult1(pHS,diffHS,alpha/2,adj)
    rejCC <- mult1(pCC,diffCC,alpha/2,adj)
    rejBach <- mult1(pBach,diffBach,alpha/2,adj)

    list(trend=c(rejHS$trend,rejBach$trend)[names(pTrendL)],
         diff=c(rejHS$diff,rejBach$diff)[names(pDiffL)])
}

mult <- function(pTrendL,pDiffL,alpha,adj){

    colnum <- ifelse(adj,2,1)
    pHS <- pTrendL$hsTot[colnum]
    pBach <- pTrendL$bachTot[colnum]

    pTrend <- do.call('rbind',pTrendL)[,colnum]
    pDiff <- do.call('rbind',pDiffL)[,colnum]

    ngroup <- length(pTrend)
    rejDiff <- rejTrend <- rep(FALSE,ngroup)

    pRank <- rank(pTrend)

    for(i in ngroup:1){
        astar <- alpha/(ngroup-i+1)
        if(pTrend[pRank==i]<= astar){
            rejTrend[pRank<=i] <- TRUE
            rejDiff[pRank<= i & pDiff<= astar] <- TRUE
            stop
        }
    }
    list(trend=toList(pTrendL,rejTrend),diff=toList(pTrendL,rejDiff))
}

mult3 <- function(pTrend,alpha,adj=TRUE){
    ## control FWER for overall trends
    ## control FDR for subgroups
    ppp <- do.call('rbind',pTrend)
    ps <- if(adj) ppp[,'adj'] else ppp[,'unadj']
    padj <- c(p.adjust(ps[1:3],method='holm'),
              p.adjust(ps[-c(1:3)],method='fdr'))
    rej <- padj<alpha
    list(trend=toList(pTrend,rej))
}

mult4 <- function(pTrend,pDiff,alpha,adj=TRUE){
    ## control FWER for overall trends
    ## control FDR for subgroup trends and hearing/deaf differences
    ppp <- do.call('rbind',pTrend)
    ps <- if(adj) ppp[,'adj'] else ppp[,'unadj']
    nTrend <- length(ps)

    ddd <- do.call('rbind',pDiff)
    psd <- if(adj) ddd[,'adj'] else ddd[,'unadj']
    ## ps <- c(ps,psd)

    padj <- c(p.adjust(ps[1:3],method='holm'),
              p.adjust(ps[-c(1:3)],method='fdr'))
  rej <- padj<alpha
  rej <- c(rej,rep(FALSE,length(psd)))
    list(trend=toList(pTrend,rej[1:nTrend]),diff=toList(pDiff,rej[(nTrend+1):length(rej)]))
}

   ##  rej,trend,nn,cc=1){
##   for(rr in rej) if(!is.finite(rr[[trend]][[nn]][cc])) return('')
##   ifelse(rej$`0.001`[[trend]][[nn]][cc],'***',
##            ifelse(rej$`0.01`[[trend]][[nn]][cc],'**',
##            ifelse(rej$`0.05`[[trend]][[nn]][cc],'*',
##                   ifelse(rej$`0.1`[[trend]][[nn]][cc],'.',''))))
## }

getSampleSizes <- function(nnn,ot,deaf){
    ccc <- combineDat(nnn,ot,deaf)
    out <- with(ccc,sapply(unique(tdat[[subCols[1]]]), function(x) round(mean(tdat$Freq[tdat[[subCols[1]]]==x]))))
    if(is.null(names(out))) names(out) <- with(ccc,unique(tdat[[subCols[1]]]))
    out
}

#### functions for plotting gaps over time
gapDat <- function(nnn,ot){
  ccc <- combineDat(nnn,ot)

  gdat <- ccc$tdat%>%
    select(-Freq)%>%
    group_by_at(vars(-DEAR,-y,-se))%>%
    summarize(gap=y[DEAR==2]-y[DEAR==1],se=sqrt(se[DEAR==1]^2+se[DEAR==2]^2))%>%
    mutate(year=year+2007)

  list(gdat=gdat,subCols=ccc$subCols,subsets=ccc$subsets)
}

plotGap <- function(gdat,errbar=TRUE,se=TRUE){

  p <- if(length(gdat$subCols)==1){
    ggplot(gdat$gdat,aes(year,gap))
  } else{
    if(gdat$subCols[2]=='SEX')
      gdat$gdat <- mutate(gdat$gdat,year=year+ifelse(SEX=='Male',1,-1)*0.1)
    ggplot(gdat$gdat,aes_string("year","gap",color=gdat$subCols[2]))
  }

  p <- p+
    geom_point()+
    geom_smooth(method='lm',se=se)+
    scale_x_continuous('Year',breaks=2008:2018,minor_breaks=NULL)+
    ylab('Hearing-Deaf Gap (Percentage Points)')

  if(length(gdat$subCols)>1) p <- p+scale_color_manual(name=NULL,values=subwayPalette)

  if(errbar) p <- p+geom_errorbar(aes(ymin=gap-2*se,ymax=gap+2*se),width=0)

  p
}

gapAdjOne <- function(tab,...){
  names(tab)[grep('edLev',names(tab))] <- 'y'
  mod <- lm_robust(y~DEAR+as.factor(AGEP),data=tab,subset=!is.na(y)&!is.na(se))
  out <- as.data.frame(summary(mod)$coef)['DEAR',1:2]
  if(any(is.nan(out[['Std. Error']]))){
    mod <- update(mod,se_type='HC1')
    out <- as.data.frame(summary(mod)$coef)['DEAR',1:2]
  }
  out
}

gapAdj <- function(dat){
  names(dat)[grep('edLev',names(dat))] <- 'y'
  subCols <- setdiff(names(dat),c('y','AGEP','se','Freq','se1','se2','year','DEAR'))

  if('SEX'%in%names(dat)) dat$SEX <- c('Male','Female')[dat$SEX]

  if(!length(subCols))
    return(gapAdjOne(dat))

  dat%>%
    group_by_at(vars(!! sym(subCols)))%>%
    group_modify(gapAdjOne)%>%
    mutate(year=dat$year[1])
}

gapAdjYr <- function(nnn,overTimeAge,weightDat=NULL){
  gaps <- lapply(overTimeAge, function(ota) gapAdj(ota[[nnn]]))

  ests <- do.call('rbind',lapply(gaps,function(x) x$Estimate))
  ses <- do.call('rbind',lapply(gaps,function(x) x$`Std. Error`))
  rownames(ests) <- rownames(ses) <- 2007+seq(length(gaps))
  if(ncol(ests)>1) colnames(ests) <- colnames(ses) <- gaps[[1]][[1]]

  gtrends <- gapTrend(nnn,overTimeAge,weightDat=weightDat)
  ests <- rbind(ests,trend=gtrends$Estimate)
  ses <- rbind(ses,trend=gtrends$`Std. Error`)
  trendPs <- gtrends$`Pr(>|t|)`
  if(ncol(ests)>1) names(trendPs) <- unlist(gtrends[,1])

  list(ests=ests,ses=ses,pvals=trendPs)
}

gapDatAdj <- function(nnn,overTimeAge){
  gay <- gapAdjYr(nnn,overTimeAge)
  ccc <- combineDat(nnn,overTimeAge)
  if(ncol(gay$ests)==1)
    return(
      list(
        gdat=data.frame(year=seq(nrow(gay$est)-1)+2007,gap=gay$ests[-nrow(gay$ests)],se=gay$se[-nrow(gay$ests)]),
        subCols=ccc$subCols,
        subsets=ccc$subsets
      )
    )

  gdat <- distinct(ccc$tdat[,c(ccc$subCols[2],'year')])
  gdat$year <- gdat$year+2007
  gdat$gap <- gdat$se <- NA
  for(i in 1:nrow(gdat)){
    gdat$gap[i] <- gay$ests[as.character(gdat$year[i]),gdat[[ccc$subCols[2]]][i]]
    gdat$se[i] <- gay$ses[as.character(gdat$year[i]),gdat[[ccc$subCols[2]]][i]]
  }

  list(
    gdat=gdat,
    subCols=ccc$subCols,
    subsets=ccc$subsets
  )
}

makeGapTab <- function(gap1,moe=TRUE){
  ests <- gap1$ests
  ses <- gap1$ses
  strs <- cbind(gap1$stars,c('',''))

  tab <- ests
  for(j in 1:ncol(ests)){
    for(i in 1:nrow(ests)){
      if(moe) tab[i,j] <- paste(round(ests[i,j],1),'$\\pm$',round(ses[i,j]*1.96,1))
      else tab[i,j] <- paste0(round(ests[i,j],1),' (',round(ses[i,j],1),')')
    }
    tab['trend',j] <- paste0(tab['trend',j],strs[j])
  }
  tab
}


text:
insert_game = function() {
    #Connect to the database
    drv <- dbDriver("SQLite")
    con = dbConnect(drv, dbfile)

    tables = dbListTables(con)

    t = 'Wisconsin'
    tt = list(team_name=t)

    if (!("teams" %in% tables)) {
        dbWriteTable(con, 'teams', as.data.frame(tt), row.names=FALSE)
    }

    dbWriteTable(con, t$name, t$dataframe)
    team_row = dbGetQuery(con, paste("SELECT rowid, team_name FROM teams WHERE team_name=='", t, "'",  sep=""))

    if (nrow(teamrow)==0) {
        dbWriteTable(con, 'teams', as.data.frame(tt), row.names=FALSE)
        team_row = dbGetQuery(con, "select last_insert_rowid()")
    }

    team_id = team_row[1]

    #Add the game to the games table
    dbWriteTable(con, "games", as.data.frame(gamemeta), row.names=FALSE)
    game_id = dbGetQuery(con, "select last_insert_rowid()")

    for (i in 1:length(drives)) {
        drives[[i]]$pbp <- NULL
        drives[[i]]$game_id = game_id
    }
}

ext: r

max_stars_repo_name: jericson/pbp

max_issues_repo_name: jericson/pbp

max_forks_repo_name: jericson/pbp


text:
#' @title Graph Information
#'
#' @description Gets graph information
#'
#' @param edgelist basically a two-column matrix with edges
#' (see \code{\link{graph}})
#' @param vertices optional vector of vertex names corresponding with
#' those in the edgelist
#' @param sorted logical to indicate if nodes should be sorted
#' (default \code{FALSE})
#' @param decreasing logical to indicate type of sorting
#' (used only when \code{sorted=TRUE})
#' @param ordering optional numeric or string vector providing the
#' ordering of nodes. When provided, this parameter overrides
#' \code{sorted=TRUE}). See the details section for more information.
#' @param labels optional string vector with labels for the nodes
#' @export
#' @keywords internal
graph_info <-
  function(edgelist, vertices, sorted = FALSE, decreasing = FALSE,
           ordering = NULL, labels = NULL)
  {
    # ======================================================
    # Checking arguments
    # ======================================================
    # edgelist as a two-column matrix
    if (!is.matrix(edgelist) || ncol(edgelist) != 2)
      stop("\nSorry, 'edgelist' must be a two column matrix")

    num_edges = nrow(edgelist)
    # get nodes (this could be numeric or character)
    if(hasArg(vertices)){
      #to deal with singleton nodes
      nodes = vertices
    }else{
      nodes = unique(as.vector(t(edgelist)))
    }
    num_nodes = length(nodes)
    # check labels (i.e. node names)
    if (!is.null(labels))
    {
      if (length(labels) != num_nodes)
        stop("\nLength of 'labels' differs from number of nodes")
    } else {
      labels = nodes
    }

    # auxiliar order (this may change if sorted or ordering required)
    aux_ord = 1:num_nodes

    # If sorted is required, ennumerate nodes
    if (sorted) {
      ordered_nodes = order(nodes, decreasing = decreasing)
      nodes = nodes[ordered_nodes]
      labels = labels[ordered_nodes]
      # auxiliar order
      aux_ord = ordered_nodes
    }

    # If ordering is provided, re-ennumerate nodes
    if (!is.null(ordering))
    {
      if (length(ordering) != num_nodes) {
        stop("\nLength of 'ordering' differs from number of nodes")
      }

      if (is.character(ordering)) {
        # make sure labels contains elements in ordering
        unmatched_ordering <- !(ordering %in% labels)
        if (any(unmatched_ordering)) {
          undetected = ordering[unmatched_ordering]
          stop(sprintf("\nUnrecognized values in ordering: '%s'", undetected))
        }
        ordering = match(ordering, labels)
      }

      nodes = nodes[ordering]
      labels = labels[ordering]
      # auxiliar order
      aux_ord = ordering
    }

    ## output
    list(
      nodes = nodes,
      labels = labels,
      num_nodes = num_nodes,
      num_edges = num_edges,
      aux_ord = aux_ord
    )
  }

#' @title X or Y coordinates of node locations
#'
#' @description
#' Gives axis locations of each node
#'
#' @param num_nodes number of nodes
#' @param aux_ord vector with the index number for ordering the nodes
#' @param labels optional string vector with labels for the nodes
#' @export
xynodes <- function(num_nodes, aux_ord, labels)
{
  # ======================================================
  # Coordinates of nodes (i.e. vertices)
  # ======================================================
  # node labels at equal distances from each other
  nf = rep(1 / num_nodes, num_nodes)
  # center coordinates of node labels
  fin = cumsum(nf)
  ini = c(0, cumsum(nf)[-num_nodes])
  centers = (ini + fin) / 2
  names(centers) = labels[aux_ord]

  # output
  centers
}

#' @title Arc Radius Locations
#'
#' @description Computes the location and radius of each arc
#'
#' @param edgelist 2-column matrix
#' @param nodes vector of nodes
#' @param centers vector with xy-positions of nodes
#' @return a list with locations and radios
#' @return \item{locs}{locations}
#' @return \item{radios}{radius values}
#' @export
#' @keywords internal
arc_radius_locs <- function(edgelist, nodes, centers)
{
  # ======================================================
  # Coordinates of arcs (i.e. edges)
  # ======================================================
  # handy matrix with numeric indices '1:FROM' , '2:TO'
  edges_from_to = matrix(0, nrow(edgelist), 2)
  for (i in 1L:nrow(edgelist))
  {
    edges_from_to[i,1] = centers[which(nodes == edgelist[i,1])]
    edges_from_to[i,2] = centers[which(nodes == edgelist[i,2])]
  }

  # maximum radius of arcs
  radios = abs(edges_from_to[,1] - edges_from_to[,2]) / 2
  max_radios = which(radios == max(radios))
  max_rad = unique(radios[max_radios] / 2)

  # arc locations
  locs = rowSums(edges_from_to) / 2

  # output
  list(locs = locs, radios = radios)
}

#' @title Above or Below
#'
#' @description Determines how arcs should be displayed
#' @details
#' If \code{horizontal = TRUE} then arcs can be plotted above or
#' below the horizontal axis \cr
#' If \code{horizontal = FALSE} then arcs can be plotted to the right or
#' left of the vertical axis
#' @param edgelist two-column matrix
#' @param above optional numeric or logical vector indicating what edges
#' (arcs) should be plotted above (or to the right of) of chosen axis
#' If \code{above = NULL} then all arcs are plotted above (or to the right)
#' If \code{above} is numeric, it cannot contain both positive and negative
#' indices.
#' If \code{above} is logical, its length must equal the number of rows in
#' \code{edgelist}
#' @return a logical vector indicating how arcs should be displayed
#' @export
#' @keywords internal
above_below <- function(edgelist, above)
{
  # ======================================================
  # Coordinates of arcs (i.e. edges) below the axis
  # ======================================================
  # check above
  if (is.null(above)) {
    above = rep(TRUE, nrow(edgelist))
  } else {
    if (length(above) > nrow(edgelist))
      stop("\nlength of 'above' exceeds number of rows in 'edgelist'")
    # check numeric above and convert to logical
    if (is.numeric(above)) {
      above_positive <- any(above > 0)
      above_negative <- any(above < 0)
      if (above_positive & above_negative)
        stop("\n'above' cannot contain both negative and positive indices")
      # convert to logical
      if (all(above > 0)) {
        above = 1:nrow(edgelist) %in% above
      }
      if (all(above < 0)) {
        above <- !(-(1:nrow(edgelist)) %in% above)
      }
      if (all(above == 0)) {
        above = rep(FALSE, nrow(edgelist))
      }
    }
    # check logical above
    if (is.logical(above)) {
      if (length(above) != nrow(edgelist))
        stop("\nlength of 'above' must equal number of rows in 'edgelist'")
    }
  }

  # output
  above
}

#' @title Minimum and Maximum Margin Limits
#' @description Computes the minimum and maximum margin limits of
#' plotting region
#' @param radios vector of arc radius
#' @param above logical vectors indicating whether arcs should be displayed
#' @return list with minimum and maximum margin limits
#' @export
#' @keywords internal
min_max_margin <- function(radios, above)
{
  # determine maximum radius
  max_radios = which(radios == max(radios))

  # minimum and maximum margin limits
  lim_min = 0
  lim_max = 0

  above_radios = radios[above]
  if (length(above_radios > 0)) {
    max_above_radios = which(above_radios == max(above_radios))[1]
    lim_max = above_radios[max_above_radios]
  }

  below_radios = radios[!above]
  if (length(below_radios > 0)) {
    max_below_radios = which(below_radios == max(below_radios))[1]
    lim_min = -1 * below_radios[max_below_radios]
  }

  # margin limits
  list(min = lim_min, max = lim_max)
}

#' @title Arc Diagram Plot
#'
#' @description
#' Give me an edgelist and I'll help you plot a pretty damn arc diagram
#'
#' @details
#' The arcs are scaled such that they fit in a plot region with its
#' x-axis ranging from zero to one. Node symbols and labels can be
#' optionally displayed. Node symbols are displayed through
#' the function \code{points}. In turn, node labels are displayed
#' through the function \code{mtext}.
#'
#' When \code{ordering} is provided in numeric format and node labels are
#' strings, the labels are alphabetically ordered first, and then nodes are
#' sorted according to the provided \code{ordering}.
#'
#' If \code{ordering} is provided in string format, the node labels must be
#' strings as well. The nodes will be sorted according to \code{ordering}.
#'
#' @param edgelist basically a two-column matrix with edges
#' (see \code{\link{graph}})
#' @param vertices optional vector of vertex names corresponding with
#' those in the edgelist
#' @param sorted logical to indicate if nodes should be sorted
#' (default \code{FALSE})
#' @param decreasing logical to indicate type of sorting
#' (used only when \code{sorted=TRUE})
#' @param ordering optional numeric or string vector providing the
#' ordering of nodes. When provided, this parameter overrides
#' \code{sorted=TRUE}). See the details section for more information.
#' @param labels optional string vector with labels for the nodes
#' @param horizontal logical indicating whether to plot
#' in horizontal orientation
#' @param above optional vector indicating which arcs should be displayed
#' above (or to the right) and below (or to the left) of the axis
#' @param col.arcs color for the arcs (default \code{"gray50"})
#' @param lwd.arcs line width for the arcs (default 1)
#' @param lty.arcs line type for the arcs (see \code{\link{par}})
#' @param lend the line end style for the arcs (see \code{\link{par}})
#' @param ljoin the line join style for the arcs (see \code{\link{par}})
#' @param lmitre the line mitre limit for the arcs (see \code{\link{par}})
#' @param show.nodes logical indicating whether to show node symbols
#' @param pch.nodes plotting 'character', i.e. symbol to use when
#' plotting nodes (\code{pch.nodes=0:25})
#' @param cex.nodes expansion of the node symbols (default 1)
#' @param col.nodes color of the node symbols (default \code{"gray50"})
#' @param bg.nodes background (fill) color for the node symbols
#' given by \code{pch.nodes=21:25}
#' @param lwd.nodes line width for drawing node symbols
#' (see \code{\link{points}})
#' @param show.labels logical indicating whether to show node labels
#' @param col.labels color of the node labels (default \code{"gray50"})
#' @param cex.labels expansion of node labels (default \code{"gray50"})
#' @param las numeric in {0,1,2,3}; the style of axis labels
#' (see \code{\link{par}})
#' @param font font used for node labels (see \code{\link{par}})
#' @param line on which margin line the node labels are displayed,
#' starting at 0 counting outwards (see \code{\link{mtext}})
#' @param outer use outer margins, if available, to plot node labels
#' (see \code{\link{mtext}})
#' @param adj adjustment for each string in reading direction
#' (see \code{\link{mtext}})
#' @param padj adjustment for each string perpendicular to
#' the reading direction (see \code{\link{mtext}})
#' @param axes logical indicating whether to plot the axes
#' (default \code{FALSE})
#' @param ... further graphical parameters (see \code{\link{par}}), including
#' \code{family}, \code{xpd}, \code{main}, \code{asp}, etc.
#' @author Gaston Sanchez
#' @seealso \code{\link{xynodes}}
#' @export
#' @examples
#'
#'  \dontrun{
#'  # create an edgelist
#'  un_graphe <- rbind(
#'  c("fromage", "pain"),
#'  c("pain", "vin"),
#'  c("vin", "biere"),
#'  c("cidre", "biere"),
#'  c("foie", "fromage"),
#'  c("pain", "foie"))
#'
#'  # deafult arcplot
#'  arcplot(un_graphe)
#'  # vertical display
#'  arcplot(un_graphe, horizontal=FALSE)
#'  # arcplot with arcs above and below axis
#'  arcplot(un_graphe, above = c(1, 3, 5))
#'  # nodes sorted alphabetically (increasing)
#'  arcplot(un_graphe, sorted=TRUE)
#'  # nodes sorted alphabetically (decreasing)
#'  arcplot(un_graphe, sorted=TRUE, decreasing = TRUE)
#'  # provided order for nodes
#'  new_order = c("vin", "biere", "cidre", "fromage", "foie", "pain")
#'  arcplot(un_graphe, ordering = new_order)
#'
#'
#'  # generate graphs
#'  ring_graph = graph.ring(10)
#'  star_graph = graph.star(10, mode="out")
#'  tree_graph = graph.tree(10, 2)
#'
#'  # add names to nodes
#'  V(ring_graph)$name = letters[1:vcount(ring_graph)]
#'  V(star_graph)$name = paste("Node", 1:vcount(star_graph))
#'  V(tree_graph)$name = paste("V", 1:vcount(tree_graph), sep='')
#'
#'  # extract edgelist
#'  ring_edges = get.edgelist(ring_graph)
#'  star_edges = get.edgelist(star_graph)
#'  tree_edges = get.edgelist(tree_graph)
#'
#'  # arc diagram
#'  arcplot(ring_edges, labels=V(ring_graph)$name, las=1)
#'  arcplot(star_edges, labels=V(star_graph)$name, las=2)
#'  arcplot(tree_edges, labels=V(tree_graph)$name, las=2)
#'
#'  # compare to plot.igraph
#'  plot(ring_graph, vertex.label=V(ring_graph)$name)
#'  plot(star_graph, vertex.label=V(star_graph)$name)
#'  plot(tree_graph, vertex.label=V(tree_graph)$name)
#'  }
#'
arcplot <- function(
  edgelist, vertices, sorted = FALSE, decreasing = FALSE, ordering = NULL,
  labels = NULL, horizontal = TRUE, above = NULL,
  col.arcs = "#5998ff77", lwd.arcs = 1.8, lty.arcs = 1,
  lend = 1, ljoin = 2, lmitre = 1, show.nodes = TRUE, pch.nodes = 19,
  cex.nodes = 1, col.nodes = "gray80", bg.nodes = "gray80", lwd.nodes = 1,
  show.labels = TRUE, col.labels = "gray55",
  cex.labels = 0.9, las = 2, font = 0, line = 0,
  outer = FALSE, adj = NA, padj = NA, axes = FALSE, xlim = NULL, ylim = NULL, ...)
{
  # Get graph information
  if (hasArg(vertices)) {
    nodes_edges = graph_info(edgelist, vertices = vertices, sorted = sorted,
                             decreasing = decreasing,
                             ordering = ordering, labels = labels)
  } else {
    nodes_edges = graph_info(edgelist, sorted = sorted, decreasing = decreasing,
                             ordering = ordering, labels = labels)
  }
  nodes = nodes_edges$nodes
  num_nodes = nodes_edges$num_nodes
  num_edges = nodes_edges$num_edges
  aux_ord = nodes_edges$aux_ord
  labels = nodes_edges$labels

  # x-y node coordinates
  centers = xynodes(num_nodes, aux_ord, labels)

  # determine above or below display of arcs
  above = above_below(edgelist, above)

  # arc radius and locations
  radios_locs = arc_radius_locs(edgelist, nodes, centers)
  radios = radios_locs$radios
  locs = radios_locs$locs

  # ======================================================
  # Graphical parameters for Arcs
  # ======================================================
  # color of arcs
  if (length(col.arcs) != num_edges)
    col.arcs = rep(col.arcs, length=num_edges)
  # line width of arcs
  if (length(lwd.arcs) != num_edges)
    lwd.arcs = rep(lwd.arcs, length=num_edges)
  # line type of arcs
  if (length(lty.arcs) != num_edges)
    lty.arcs = rep(lty.arcs, length=num_edges)

  # ======================================================
  # Graphical parameters for Nodes
  # ======================================================
  # pch symbol of nodes
  if (length(pch.nodes) != num_nodes) {
    pch.nodes = rep(pch.nodes, length = num_nodes)
  }
  pch.nodes = pch.nodes[aux_ord]
  # cex of nodes
  if (length(cex.nodes) != num_nodes) {
    cex.nodes = rep(cex.nodes, length = num_nodes)
  }
  cex.nodes = cex.nodes[aux_ord]
  # color of nodes
  if (length(col.nodes) != num_nodes) {
    col.nodes = rep(col.nodes, length = num_nodes)
  }
  col.nodes = col.nodes[aux_ord]
  # bg of nodes
  if (length(bg.nodes) != num_nodes) {
    bg.nodes = rep(bg.nodes, length = num_nodes)
  }
  bg.nodes = bg.nodes[aux_ord]
  # line widths of nodes
  if (length(lwd.nodes) != num_nodes) {
    lwd.nodes = rep(lwd.nodes, length = num_nodes)
  }
  lwd.nodes = lwd.nodes[aux_ord]

  # ======================================================
  # Graphical parameters for Node Labels
  # ======================================================
  # color of labels
  if (length(col.labels) != num_nodes) {
    col.labels = rep(col.labels, length = num_nodes)
  }
  col.labels = col.labels[aux_ord]
  # cex of labels
  if (length(cex.labels) != num_nodes) {
    cex.labels = rep(cex.labels, length = num_nodes)
  }
  cex.labels = cex.labels[aux_ord]

  # ======================================================
  # Plot arc diagram (horizontally or vertically)
  # ======================================================
  # auxiliar vector for plotting arcs
  z = seq(0, pi, length.out = 100)

  if (horizontal) {
    side = 1
  } else {
    side = 2
  }

  if (is.null(xlim)) {
    if (horizontal) {
      xlim = c(-0.015, 1.015)
      x_nodes = centers
    } else {
      xlims = min_max_margin(radios, above)
      xlim = c(xlims$min, xlims$max)
      x_nodes = rep(0, num_nodes)
    }
  } else {
    if (horizontal) {
      x_nodes = centers
    } else {
      x_nodes = rep(0, num_nodes)
    }
  }

  if (is.null(ylim)) {
    if (horizontal) {
      ylims = min_max_margin(radios, above)
      ylim = c(ylims$min, ylims$max)
      y_nodes = rep(0, num_nodes)
    } else {
      ylim = c(-0.015, 1.015)
      y_nodes = centers
    }
  } else {
    if (horizontal) {
      y_nodes = rep(0, num_nodes)
    } else {
      y_nodes = centers
    }
  }

  # open empty plot window
  plot(0.5, 0.5, xlim = xlim, ylim = ylim, type = "n",
       xlab = "", ylab = "", axes = axes, ...)
  # add each edge
  for (i in 1L:num_edges)
  {
    # get radius length
    radio = radios[i]
    if (horizontal) {
      # x-y coords of each arc
      x_arc = locs[i] + radio * cos(z)
      if (above[i]) { # above axis
        y_arc = radio * sin(z)
      } else {  # below axis
        y_arc = radio * sin(-z)
      }
    } else {
      # x-y coords of each arc
      y_arc = locs[i] + radio * cos(z)
      if (above[i]) { # above axis
        x_arc = radio * sin(z)
      } else {  # below axis
        x_arc = radio * sin(-z)
      }
    }

    # plot arc connecting nodes
    lines(x_arc, y_arc, col=col.arcs[i], lwd=lwd.arcs[i], lty=lty.arcs[i],
          lend=lend, ljoin=ljoin, lmitre=lmitre)
    # add node symbols with points
    if (show.nodes) {
      points(x=x_nodes, y=y_nodes, pch=pch.nodes,
             col=col.nodes, bg=bg.nodes, cex=cex.nodes, lwd=lwd.nodes)
    }
    # add node labels with mtext
    if (show.labels) {
      mtext(labels, side=side, line=line, at=centers, cex=cex.labels, outer=outer,
            col=col.labels, las=las, font=font, adj=adj, padj=padj, ...)
    }
  }
}

#' @title Node Coordinates
#'
#' @description
#' Computes axis locations of each node. This function can be helpful when
#' you want to separately plot the node labels using the function mtext.
#'
#' @param edgelist basically a two-column matrix with edges
#' (see \code{\link{graph}})
#' @param sorted logical to indicate if nodes should be sorted
#' @param decreasing logical to indicate type of sorting
#' (used only when \code{sorted=TRUE})
#' @param ordering optional numeric vector providing the ordering of nodes
#' (when provided, this parameter overrides \code{sorted=TRUE})
#' @param labels character vector with labels for the nodes
#' @return a vector with the location of nodes in the x-axis
#' @author Gaston Sanchez
#' @seealso \code{\link{arcplot}}
#' @export
#' @examples
#'
#'  \dontrun{
#'  # generate a graph
#'  some_graph = graph.ring(10)
#'
#'  # add names to nodes
#'  V(some_graph)$name = letters[1:vcount(some_graph)]
#'
#'  # extract edgelist
#'  edgelist = get.edgelist(some_graph)
#'
#'  # (default) arc diagram
#'  arcplot(edgelist, labels=V(some_graph)$name, las=1)
#'
#'  # get x-axis coordinates of nodes
#'  xcoords = node_coords(edgelist, labels=V(some_graph)$name)
#'
#'  # arc diagram with various labels
#'  arcplot(edgelist, show.labels=FALSE, show.nodes=TRUE)
#'  mtext(V(some_graph)$name, side=1, line=0, at=xcoords)
#'  mtext(rep("node",10), side=1, line=1, at=xcoords, col="gray90")
#'  }
#'
node_coords <- function(
  edgelist, sorted = FALSE, decreasing = FALSE, ordering = NULL,
  labels = NULL)
{
  # Get graph information
  nodes_edges = graph_info(edgelist, sorted = sorted, decreasing = decreasing,
                           ordering = ordering, labels = labels)

  nodes = nodes_edges$nodes
  num_nodes = nodes_edges$num_nodes
  num_edges = nodes_edges$num_edges
  aux_ord = nodes_edges$aux_ord
  labels = nodes_edges$labels

  # x-y node coordinates
  centers = xynodes(num_nodes, aux_ord, labels)
}


text:
library(shiny)
library(plyr)
#library(tidyr)
library(dplyr)
require(devtools)
#install_github('ramnathv/rCharts')
require(rCharts)

# Define server logic for slider examples
shinyServer(function(input, output,session) {

  DatasetInput<-reactive({
    inFile<-input$file1
    if (is.null(inFile))
      return(NULL)

    read.csv(inFile$datapath, header=input$header, sep=input$sep,
             quote=input$quote)
  })

  #first tab test
  output$summary1<-renderPrint({
    dataset<-DatasetInput()
    if(is.null(dataset)) return ("For Infection,Movement,Stifness Summary by Windows")
    #summary(dataset)
    a<-plyr::count(dataset,"Window", "Infection")
    a<-plyr::rename(a,replace=c("freq"="Infection"))
    b<-plyr::count(dataset,"Window", "Movement")
    b<-plyr::rename(b,replace=c("freq"="Movement"))
    c<-plyr::count(dataset,"Window", "Stiffness")
    c<-plyr::rename(c,replace=c("freq"="Stiffness"))
    d<-plyr::count(dataset,"Window")
    # d<-plyr::rename(d,replace=c("x"="Window"))
    ab<-plyr::join(x=a,y=b,by="Window",type="left")
    abc<-plyr::join(x=ab,y=c,by="Window",type="left")
    abcd<-plyr::join(x=abc,y=d,by="Window",type="left")
    abcd
  })

  values<-reactiveValues(data = NULL)
  observeEvent(input$file1, {
    values$data <- read.csv(input$file1$datapath, header=input$header, sep=input$sep,
                            quote=input$quote)
    write.csv(values$data,'updated.csv') #write out file will be under Rshiny folder
  })

  vals<-eventReactive(input$save,{

    #summary(dataset)
    if (input$x=="Yes") {xx=1} else {xx=0}
    if (input$y=="Yes") {yy=1} else {yy=0}
    if (input$z=="Yes") {zz=1} else {zz=0}

    values$data[nrow(values$data)+1,1]<-input$ID
    values$data[nrow(values$data),2]<-xx
    values$data[nrow(values$data),3]<-yy
    values$data[nrow(values$data),4]<-zz
    values$data[nrow(values$data),5]<-input$days
    temp <- values$data
    values$data <- temp
    write.csv(values$data,'updated.csv') #write out file will be under Rshiny folder
    values$data[,1:5]

  })

  output$summary2<-renderPrint({

    #vals()

    dataset<-vals()
    #if(is.null(dataset)) return ("For Infection,Movement,Stifness Summary by Windows")
    #summary(dataset)
    a<-plyr::count(dataset,"Window", "Infection")
    a<-plyr::rename(a,replace=c("freq"="Infection"))
    b<-plyr::count(dataset,"Window", "Movement")
    b<-plyr::rename(b,replace=c("freq"="Movement"))
    c<-plyr::count(dataset,"Window", "Stiffness")
    c<-plyr::rename(c,replace=c("freq"="Stiffness"))
    d<-plyr::count(dataset,"Window")
    # d<-plyr::rename(d,replace=c("x"="Window"))
    ab<-plyr::join(x=a,y=b,by="Window",type="left")
    abc<-plyr::join(x=ab,y=c,by="Window",type="left")
    abcd<-plyr::join(x=abc,y=d,by="Window",type="left")
    abcd

  })

  output$table1<-renderTable({
    head(DatasetInput(),n=input$obs)
  })

  # Reactive expression to compose a data frame containing all of
  # the values
  TableValues2 <- reactive({

    # Compose data frame
    data.frame(
      Name = c("P",
               "Q",
               "R"
      ),
      Value = as.character(c(input$x,
                             input$y,
                             input$z)),
      stringsAsFactors=FALSE)
  })

  # Show the values using an HTML table
  #   output$values2 <- renderTable({
  #     TableValues2()
  #   })

  TableValues3 <- reactive({

    # Compose data frame
    data.frame(
      Name = c("Stochastic walk"
      ),
      Window1 = as.character( stage2()[2]),
      Window2=as.character( stage2()[1]),
      stringsAsFactors=FALSE)
  })

  #second tab contents
  output$values3 <- renderTable({
    TableValues3()
  })

  updates<-eventReactive(input$ShowPlot,{dataset<-read.csv("updated.csv", header=TRUE,  sep = ",",stringsAsFactors = FALSE)

  #if(is.null(dataset)) return ("Waiting for dataset")
  #summary(dataset)
  a<-plyr::count(dataset,"Window", "Infection")
  a<-plyr::rename(a,replace=c("freq"="Infection"))
  b<-plyr::count(dataset,"Window", "Movement")
  b<-plyr::rename(b,replace=c("freq"="Movement"))
  c<-plyr::count(dataset,"Window", "Stiffness")
  c<-plyr::rename(c,replace=c("freq"="Stiffness"))
  d<-plyr::count(dataset,"Window")
  #d<-plyr::rename(d,replace=c("x"="Window"))
  ab<-plyr::join(x=a,y=b,by="Window",type="left")
  abc<-plyr::join(x=ab,y=c,by="Window",type="left")
  abcd<-plyr::join(x=abc,y=d,by="Window",type="left")
  library(tidyr)
  abcdLong<-as.data.frame(abcd%>%tidyr::gather(status,freq1, Infection:freq))
  })

  output$myChart <- renderChart({
    p1 <- nPlot(freq1 ~ Window, group = 'status', data = updates(),
                type = input$type, dom = 'myChart', width = 800)
    #p6$chart(color = c('brown', 'blue', '#594c26', 'green'), stacked = input$stack)
    return(p1)
  })

  stage2<-eventReactive(input$stage2,{
    source("Stage2backendforMac.R", local=TRUE)
    #L1<-data.frame(c("1","2","3","4","5"),L,stringsAsFactors = FALSE)
    c(Windownext, Movestep,L)  #L has 5 value

  })

  output$summary3<-renderPrint({
    stage2()
  })

  output$myChart1 <- renderChart2({
    #L has 5 value corresponding to [3:7] in vector stage2
    L1<-data.frame(c("1","2","3","4","5"),stage2()[3:7],stringsAsFactors = FALSE)

    colnames(L1)<-c("window","value")
    m1<-mPlot(x="window",y="value", data=L1,ymax=round(max(as.numeric(L1$value)),3),ymin=round(min(as.numeric(L1$value)),3),type="Line",parseTime=FALSE)
    m1$set( labels = "value" )
    return(m1)
  })

  #x    <- faithful[, 2]  # Old Faithful Geyser data
  #bins <- seq(min(x), max(x), length.out = input$bins + 1)

  # draw the histogram with the specified number of bins
  #hist(x, breaks = bins, col = 'darkgray', border = 'white')
})

max_stars_repo_name: czqiu/YiyiClinicalTrial

max_issues_repo_name: czqiu/YiyiClinicalTrial

max_forks_repo_name: czqiu/YiyiClinicalTrial


text:
# The following defines a function which has a function inside its body. This is
# *different* from a nested function; that is, `f` is different from
#
#   f1 = function () { identity }
#
# … because `f1` contains a *name* in its body. Likewise, it is different from
#
#   f2 = function () { function (x) x }
#
# … because `f2` contains a *call expression* which, when executed, *defines* a
# function, in its body. That is:
#
#   class(body(f)[[2L]])  == 'function'
#   class(body(f1)[[2L]]) == 'name'
#   class(body(f2)[[2L]]) == 'call'
f = function () { NULL }
body(f)[[2L]] = identity

#' `g` clearly isn’t a generic even though `UseMethod` is used inside a nested
#' function in its body.
#' @export
g = function (x) {
    nested = function () UseMethod('nested')
    nested()
}

#' … nor is `h`.
#' @export
h = function () { NULL }
body(h)[[2L]] = function () UseMethod('foo')

max_stars_repo_name: klmr/box

max_issues_repo_name: klmr/box

max_forks_repo_name: klmr/box


text:
###Data pre-processing###

options(warn=-1) #surpress warnings

#imports for data wrangling
library(data.table)
library(dtplyr)
library(dplyr)

#get the data - nrows set to 10000 to keep runtime manageable.
#one expansion option would be to select a time frame to visualized
NumericData <- fread("../input/train_numeric.csv", select = c("Id", "Response"),nrows = 10000)
Date <- fread("../input/train_date.csv", nrows = 10000)

#for each job identify which stations are passed through and for those store the minimum time
for (station in paste0("S",0:51))
{
  cols = min(which((grepl(station,colnames(Date)))))
  if(!cols==Inf){
    Date[,paste0(station) := Date[,cols,with = FALSE]]
  }
}

#limit data to only when passed through station X
Stations = Date[,!grepl("L",colnames(Date)),with=F]

#melt data to go from wide to long format
Stations_long = melt(Stations,id.vars=c("Id"))

#join with numeric to have Response
Stations_long <- left_join(NumericData, by="Id", Stations_long)

#remove NA entries - these are plentiful as after melting each station-job combination has its own row
Stations_filtered<-filter(Stations_long,!is.na(value))

#sort entries by ascending time
Stations_filtered<-arrange(Stations_filtered, value)

## Help sought from: https://www.kaggle.com/gingerman/bosch-production-line-performance/shopfloor-visualization

#imports for plotting
require(GGally)
library(network)
library(sna)
library(ggplot2)

#plotting format
options(repr.plot.width=5, repr.plot.height=15)

#for each row obtain the subsequent station
edgelistsComplete<-mutate(group_by(Stations_filtered, Id), nextStation = lead(variable))

#for each id find the first node to be entered
startingPoints<-select(ungroup(filter(group_by(edgelistsComplete,Id),!(variable %in% nextStation))),variable, Response)

#prior to each starting point insert an edge from a common origin
colnames(startingPoints) = c("nextStation","Response")
startingPoints$variable = "S"
paths<-select(edgelistsComplete,variable,nextStation,Response)

#for each id find the row where there is no next station (last station to be visited)
#fill this station with Response value
paths[is.na(nextStation)]$nextStation = paste("Result",paths[is.na(nextStation)]$Response)

#combine data
paths = rbind(startingPoints,paths)
paths = select(paths,-Response)
paths$nextStation = as.character(paths$nextStation)
paths$variable = as.character(paths$variable)

#rename columns for plotting
colnames(paths) <- c("Target","Source")

#flip columns in a costly way because ggnet is a little dumb and I am lazy
pathshelp = select(paths,Source)
pathshelp$Target = paths$Target
paths=pathshelp

#create network from edgelist
net = network(as.data.frame(na.omit(paths)),
              directed = TRUE)

#create a station-line mapping lookup
LineStations = NULL
for (station in unique(paths$Source)){
  if(station!="S")
  {
  x=paste0("_",station,"_")
  y=head(colnames(Date)[which(grepl(x,colnames(Date)))],1)
  y=strsplit(y,"_")[[1]][1]
  LineStations = rbind(LineStations,data.frame(Node=station,Line=y))
  }
}
LineStations = rbind(LineStations,data.frame(Node=c("Result 1","Result 0","S"),Line=c("Outcome","Outcome","START")))

#merge station-line mapping into graph for coloring purposes
x = data.frame(Node = network.vertex.names(net))
x = merge(x, LineStations, by = "Node", sort = FALSE)$Line
net %v% "line" = as.character(x)

#setup station coordinates
nodeCoordinates=data.frame(label=c("S","S0","S1","S2","S3","S4","S5","S6",
                                  "S7","S8","S9","S10","S11","S12","S13",
                                  "S14","S15","S16","S17","S18","S19",
                                  "S20","S21","S22","S23","S24","S25",
                                  "S26","S27","S28","S29","S30","S31",
                                  "S32","S33","S34","S35","S36","S37",
                                  "S38","S39","S40","S41","S43",
                                  "S44","S45","S47","S48","S49",
                                  "S50","S51","Result 0","Result 1"),
                           y=c(0,
                               1,2,3,3,4,4,5,5,6,7,7,7,
                               1,2,3,3,4,4,5,5,6,7,7,7,
                               6,6,7,7,7,
                               8,9,10,10,10,11,11,12,13,14,
                               8,9,10,11,11,12,13,14,15,15,16,
                               17,17),
                           x=c(5,
                               9,9,10,8,10,8,10,8,7,10,9,8,
                               5,5,6,4,6,4,6,4,5,6,5,4,
                               2,0,2,1,0,
                               7,7,8,7,6,8,6,7,7,7,
                               3,3,3,4,2,3,3,3,4,2,3,
                               7,3))

nodeCoordinates$y = -3 * nodeCoordinates$y

#setup initial plot
network = ggnet2(net)

#grab node list from initial plot and attach coordinates
Coordinates = select(network$data,label)
Coordinates = left_join(Coordinates,nodeCoordinates,by = "label")
Coordinates = as.matrix(select(Coordinates,x,y))

#setup plot with manual layout
network = ggnet2(net,
                 alpha = 0.75, size = "indegree",
                 label = T, size.cut = 4,
                 color = "line",palette = "Set1",
                 mode = Coordinates,
                 edge.alpha = 0.5, edge.size = 1,
                 legend.position = "bottom")

#output plot on graphics device
print(network)


text:
#' @param facet  A list of facet names used to retrieve the 100 most frequent
#' values for a field. Allowed facets are: datasetKey, highertaxonKey, rank,
#' status, extinct, habitat, and nameType. Additionally threat and
#' nomenclaturalStatus are legal values but not yet implemented, so data will
#' not yet be returned for them.
#' @param facetMincount Used in combination with the facet parameter. Set
#' \code{facetMincount={#}} to exclude facets with a count less than {#}, e.g.
#' http://bit.ly/1bMdByP only shows the type value 'ACCEPTED' because the other
#' statuses have counts less than 7,000,000
#' @param facetMultiselect  Used in combination with the facet parameter. Set
#' facetMultiselect=true to still return counts for values that are not
#' currently filtered, e.g. http://bit.ly/19YLXPO still shows all status
#' values even though status is being filtered by \code{status=ACCEPTED}

max_stars_repo_name: damianooldoni/rgbif

max_issues_repo_name: damianooldoni/rgbif

max_forks_repo_name: damianooldoni/rgbif


text:
#' ggplot2 stub
#'
#' @return ggplot2 stub
#' @export
gg <- function() {
  writeLines('ggplot(df, aes(x, y)) +\ngeom_point()')
}


text:
REBOL [
    System: "REBOL [R3] Language Interpreter and Run-time Environment"
    Title: "REBOL 3 Mezzanine: Help"
    Rights: {
        Copyright 2012 REBOL Technologies
        REBOL is a trademark of REBOL Technologies
    }
    License: {
        Licensed under the Apache License, Version 2.0
        See: http://www.apache.org/licenses/LICENSE-2.0
    }
]

spec-of: function [
    {Generate a block which could be used as a "spec block" from an action.}

    return: [block!]
    action [action!]
][
    meta: try match object! meta-of :action

    specializee: try match action! select meta 'specializee
    adaptee: try match action! select meta 'adaptee
    original-meta: try match object! any [
        meta-of :specializee
        meta-of :adaptee
    ]

    return collect [
        keep/line ensure [<opt> text!] any [
            select meta 'description
            select original-meta 'description
        ]

        return-type: try ensure [<opt> block!] any [
            select meta 'return-type
            select original-meta 'return-type
        ]
        return-note: try ensure [<opt> text!] any [
            select meta 'return-note
            select original-meta 'return-note
        ]
        if any [return-type return-note] [
            keep compose [
                return: (opt return-type) (opt return-note)
            ]
        ]

        types: try ensure [<opt> frame!] any [
            select meta 'parameter-types
            select original-meta 'parameter-types
        ]
        notes: try ensure [<opt> frame!] any [
            select meta 'parameter-notes
            select original-meta 'parameter-notes
        ]

        for-each param parameters of :action [
            keep compose [
                (param) (select types param) (select notes param)
            ]
        ]
    ]
]

description-of: function [
    {One-line summary of a value's purpose}

    return: [<opt> text!]
    v [<blank> any-value!]
][
    opt switch type of get/any 'v [
        void! [blank]
        any-array! [spaced ["array of length:" length of v]]
        image! [spaced ["size:" v/size]]
        datatype! [
            spec: ensure object! spec of v  ; "type specs" need simplifying
            copy spec/title
        ]
        action! [
            try if meta: meta-of :v [
                copy get 'meta/description
            ]
        ]
        gob! [spaced ["offset:" v/offset "size:" v/size]]
        object! [mold words of v]
        typeset! [mold make block! v]
        port! [mold reduce [v/spec/title v/spec/ref]]
        default [blank]
    ]
]

browse: function [
    "stub function for browse* in extensions/process/ext-process-init.reb"

    return: <void>
    location [<blank> url! file!]
][
    print "Browse needs redefining"
]

help: function [
    "Prints information about words and values (if no args, general help)."

    return: <void>
    'topic [<end> any-value!]
        "WORD! whose value to explain, or other HELP target (try HELP HELP)"
    /doc
        "Open web browser to related documentation."
][
    if undefined? 'topic [
        print "#[void] is a literal VOID! value"
        return
    ]

    if null? :topic [
        ;
        ; Was just `>> help` or `do [help]` or similar.
        ; Print out generic help message.
        ;
        print trim/auto copy {
            Use HELP to see built-in info:

                help insert

            To search within the system, use quotes:

                help "insert"

            To browse online topics:

                help #compiling

            To browse online documentation:

                help/doc insert

            To view words and values of a context or object:

                help lib    - the runtime library
                help self   - your user context
                help system - the system object
                help system/options - special settings

            To see all words of a specific datatype:

                help object!
                help action!
                help datatype!

            Other debug helpers:

                docs - open browser to web documentation
                dump - display a variable and its value
                probe - print a value (molded)
                source func - show source code of func
                trace - trace evaluation steps
                what - show a list of known functions
                why - explain more about last error (via web)

            Other information:

                about - see general product info
                bugs - open GitHub issues website
                changes - show changelog
                chat - open GitHub developer forum
                install - install (when applicable)
                license - show user license
                topics - open help topics website
                upgrade - check for newer versions
                usage - program cmd line options
        }
        return
    ]

    ; !!! R3-Alpha permitted "multiple inheritance" in objects, in the sense
    ; that it would blindly overwrite fields of one object with another, which
    ; wreaked havoc on the semantics of functions in unrelated objects.  It
    ; doesn't work easily with derived binding, and doesn't make a lot of
    ; sense.  But it was used here to unify the lib and user contexts to
    ; remove potential duplicates (even if not actually identical).  This
    ; does that manually, review.
    ;
    make-libuser: does [
        libuser: copy system/contexts/lib
        for-each [key val] system/contexts/user [
            if not void? get/any 'val [
               append libuser reduce [key :val]
            ]
        ]
        libuser
    ]

    switch type of :topic [
        issue! [
            ; HELP #TOPIC will browse r3n for the topic

            browse join-all [https://r3n.github.io/topics/ as text! topic]
            print newline
            return
        ]

        text! [
            ; HELP "TEXT" wildcard searches things w/"TEXT" in their name

            if types: summarize-obj/match make-libuser :topic [
                print "Found these related words:"
                for-each line sort types [
                    print line
                ]
            ] else [
                print ["No information on" topic]
            ]
            return
        ]

        path! word! [
            ; PATH! and WORD! fall through for help on what they look up to

            all [
                word? topic
                null? binding of topic
            ] then [
                print [topic "is an unbound WORD!"]
                return
            ]

            switch type of value: get/any topic [
                null [
                    print [topic "is null"]
                ]
                void! [
                    print [topic "is not defined (e.g. has a VOID! value)"]
                ]
            ] then [
                return
            ]
            enfixed: did all [action? :value | enfixed? :value]
        ]
    ] else [
        ; !!! There may be interesting meanings to apply to things like
        ; `HELP 2`, which could be distinct from `HELP (1 + 1)`, or the
        ; same, and could be distinct from `HELP :(1)` or `HELP :[1]` etc.
        ; For the moment, it just tells you the type.

        if free? :topic [
            print ["is a freed" mold type of :topic]
        ] else [
            print [mold topic "is" an mold type of :topic]
        ]
        return
    ]

    ; Open the web page for it?
    if doc and [match [action! datatype!] :value] [
        item: form :topic
        if action? get :topic [
            ;
            ; !!! The logic here repeats somewhat the same thing that is done
            ; by TO-C-NAME for generating C identifiers.  It might be worth it
            ; to standardize how symbols are textualized for C with what the
            ; documentation uses (though C must use underscores, not hyphen)
            ;
            for-each [a b] [
                "!" "-ex"
                "?" "-q"
                "*" "-mul"
                "+" "-plu"
                "/" "-div"
                "=" "-eq"
                "<" "-lt"
                ">" "-gt"
                "|" "-bar"
            ][
                replace/all item a b
            ]

            browse join-all [
                https://github.com/gchiu/reboldocs/blob/master/
                item
                %.MD
            ]
        ] else [
            remove back tail of item  ; it's a DATATYPE!, so remove the !
            browse join-all [
                http://www.rebol.com/r3/docs/datatypes/
                item
                tmp: %.html
            ]
        ]
    ]

    if datatype? :value [
        if instances: summarize-obj/match make-libuser :value [
            print ["Found these" (uppercase form topic) "words:"]
            for-each line instances [
                print line
            ]
        ] else [
            print [topic {is a datatype}]
        ]
        return
    ]

    if not action? :value [
        print collect [
            keep [(uppercase mold topic) "is" an (mold type of :value)]
            if free? :value [
                keep "that has been FREEd"
            ] else [
                keep "of value:"
                if match [object! port!] value [
                    keep newline
                    for-each line summarize-obj value [
                        keep line
                        keep newline
                    ]
                ] else [
                    keep mold value
                ]
            ]
        ]
        return
    ]

    ; The HELP mechanics for ACTION! are more complex in Ren-C due to the
    ; existence of function composition tools like SPECIALIZE, CHAIN, ADAPT,
    ; HIJACK, etc.  Rather than keep multiple copies of the help strings,
    ; the relationships are maintained in META-OF information on the ACTION!
    ; and are "dug through" in order to dynamically inherit the information.
    ;
    ; Code to do this evolved rather organically, as automatically generating
    ; help for complex function derivations is a research project in its
    ; own right.  So it tends to break--test it as much as possible.

    print "USAGE:"

    args: _  ; plain arguments
    refinements: _  ; refinements and refinement arguments

    parse parameters of :value [
        copy args any [word! | get-word! | lit-word! | issue!]
        copy refinements any path!  ; !!! Refinements may become intermixed!
        end
    ]

    ; Output exemplar calling string, e.g. LEFT + RIGHT or FOO A B C
    ; !!! Should refinement args be shown for enfixed case??
    ;
    if enfixed and [not empty? args] [
        print [_ _ _ _ args/1 (uppercase mold topic) next args]
    ] else [
        print [_ _ _ _ (uppercase mold topic) args refinements]
    ]

    meta: try meta-of :value

    print newline

    print "DESCRIPTION:"
    print [_ _ _ _ (:meta/description or '{(undocumented)})]
    print [_ _ _ _ (uppercase mold topic) {is an ACTION!}]

    print-args: function [list /indent-words] [
        for-each param list [
            type: ensure [<opt> block!] (
                select try :meta/parameter-types to-word param
            )
            note: ensure [<opt> text!] (
                select try :meta/parameter-notes to-word param
            )

            print [_ _ _ _ param (if type [unspaced ["[" type "]"]])]
            if note [
                print [_ _ _ _ _ _ _ _ note]
            ]
        ]
    ]

    print newline
    print [
        "RETURNS:"
        either :meta/return-type [mold :meta/return-type] ["(undocumented)"]
    ]
    if :meta/return-note [
        print [_ _ _ _ meta/return-note]
    ]

    if not empty? args [
        print newline
        print "ARGUMENTS:"
        print-args args
    ]

    if not empty? refinements [
        print newline
        print "REFINEMENTS:"
        print-args/indent-words refinements
    ]
]

source: function [
    "Prints the source code for an ACTION! (if available)"

    return: <void>
    'arg [word! path! action! tag!]
][
    switch type of :arg [
        tag! [
            f: copy "unknown tag"
            for-each location words of system/locale/library [
                if location: select load get location arg [
                    f: form location/1
                    break
                ]
            ]
        ]

        word! path! [
            name: arg
            f: get arg else [
                print [name "is not set to a value"]
                return
            ]
        ]
    ] else [
        name: "anonymous"
        f: :arg
    ]

    case [
        match [text! url!] :f [
            print f
        ]
        not action? :f [
            print [name "is" an mold type of :f "and not an ACTION!"]
        ]
    ] then [
        return
    ]

    ; ACTION!
    ; The system doesn't preserve the literal spec, so it must be rebuilt
    ; from combining the the META-OF information.

    write-stdout unspaced [
        mold name ":" _ "make action! [" _ mold spec-of :f
    ]

    ; While all interfaces as far as invocation is concerned has been unified
    ; under the single ACTION! interface, the issue of getting things like
    ; some kind of displayable "source" would have to depend on the dispatcher
    ; used.  For the moment, BODY OF hands back limited information.  Review.
    ;
    switch type of (body: body of :f) [
        block! [  ; FUNC, FUNCTION, PROC, PROCEDURE or (DOES of a BLOCK!)
            print [mold body "]"]
        ]

        frame! [  ; SPECIALIZE (or DOES of an ACTION!)
            print [mold body "]"]
        ]
    ] else [
        print "...native code, no source available..."
    ]
]

what: function [
    {Prints a list of known actions}

    return: [void! block!]
    'name [<end> word! lit-word!]
        "Optional module name"
    /args
        "Show arguments not titles"
    /as-block
        "Return data as block"
][
    list: make block! 400
    size: 0

    ; copy to get around error: "temporary hold for iteration"
    ctx: copy all [set? 'name try select system/modules :name ] else [lib]

    for-each [word val] ctx [
        ; word val
        if (action? :val) [
            arg: either args [
                mold parameters of :val
            ][
                description-of :val else ["(undocumented)"]
            ]
            append list reduce [word arg]
            size: max size length of to-text word
        ]
    ]
    list: sort/skip list 2

    name: make text! size
    either not as-block [
        for-each [word arg] list [
            append/dup clear name #" " size
            change name word
            print [
                name
                :arg
            ]
        ]
    ][
         list
    ]
]

pending: does [
    comment "temp function"
    print "Pending implementation."
]

bugs: func [return: <void>] [
    "View bug database."
][
    browse https://github.com/metaeducation/ren-c/issues
]

chat: func [
    "Open REBOL/ren-c developers chat forum"
    return: <void>
][
    browse http://chat.stackoverflow.com/rooms/291/rebol
]

; temporary solution to ensuring scripts run on a minimum build
;
require-commit: function [
    "checks current commit against required commit"

    return: <void>
    commit [text!]
][
    c: select system/script/header 'commit else [return]

    ; If we happen to have commit information that includes a date, then we
    ; can look at the date of the running Rebol and know that a build that is
    ; older than that won't work.
    ;
    all [
        date: select c 'date
        system/build < date

        fail [
            "This script needs a build newer or equal to" date
            "so run `upgrade`"
        ]
    ]

    ; If there's a specific ID then assume that if the current build does not
    ; have that ID then there *could* be a problem.
    ;
    if (id: select c 'id) and [id <> commit] [
        print [
            "This script has only been tested again commit" id LF

            "If it doesn't run as expected"
            "you can try seeing if this commit is still available" LF

            "by using the `do <dl-renc>` tool and look for"
            unspaced [
                "r3-" copy/part id 7 "*"
                if find-last form system/version "0.3.4" [%.exe]
            ]
        ]
    ]
]

max_stars_repo_name: giuliolunati/ren-c

max_issues_repo_name: giuliolunati/ren-c

max_forks_repo_name: giuliolunati/ren-c


text:
#' Stack custom icons with group coloring
#' @param D Data Frame
#' @param var categorical variable that is to be stacked
#' @param group grouping categorical variable that is to color the icons
#' @param icon link to any image icon on the internet. For coloring to work, has to be a PNG with transparent background
#' @param size icon size
#' @param asp icon aspect ratio
#' @param size icon size
#' @return a ggplot object that can be modified downstream
#' @export

istack = function(D, var, group = NULL, icon, size = 0.03, asp = 1, sort = TRUE) {

  # incase this is a tibble
  D = data.frame(D)

  # giving nicknames
  D['var'] = factor(D[[var]])

  if (!is.null(group)) {
    D['group'] = factor(D[[group]])
  }

  # sort the levels of variable column
  if (sort) {
    D['var'] = factor(D[['var']], names(sort(table(D[['var']]), decreasing = FALSE)))
  }

  # calculate coordinates
  D =
    do.call(
      rbind,
      lapply(
        split(D, D['var']),
        function(df) {
          if (!is.null(group)) {
            df = df[order(df['group']),]
          }
          df['n'] = seq(nrow(df))
          return(df)
        }
      )
    )

  D['image'] = icon

  # create the plot
  if (!is.null(group)) {
    base = ggplot2::ggplot(data = D, aes(x = n, y = var, color = group)) + ggplot2::labs(color = group)
  } else {
    base = ggplot2::ggplot(data = D, aes(x = n, y = var))
  }

  p = base +
    ggimage::geom_image(aes(image=image), asp = asp, size = size) +
    ggplot2::theme(
      plot.title = element_text(hjust = 0.5),
      panel.background = element_blank(),
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      axis.ticks = element_blank(),
      axis.line = element_blank()
    ) +
    ggplot2::ylab(var) +
    ggplot2::xlab('') +
    ggplot2::scale_x_continuous(breaks = scales::pretty_breaks(4), expand = c(0,0.6))

  return(p)
}

#' Simulate cancer demo dataset
#' @export
simulate_cancers = function(n = 10) {
  diseases = c(rep('Breast cancer', 3), rep('Prostate cancer', 3),
               'Brain cancer', rep('Colorectal cancer', 2), rep('Colon cancer', 3),
               'Pancreatic cancer', 'Thyroid cancer', 'Lung cancer',
               'Bladder cancer', 'Ovarian cancer', 'Sarcoma', 'Leukemia')
  treatments = c('Untreated', rep('Radiotherapy', 3), rep('Chemotherapy', 5),
                 'Targeted therapy')
  nums = sample((round(n/2)):(n*3), length(diseases), replace=T)

  cancers = data.frame(unname(do.call(c, mapply(function(d, n) {rep(d, n)}, diseases, nums))))
  colnames(cancers) = 'Disease'
  cancers['Treatment'] = sample(treatments, nrow(cancers), replace=T)
  return(cancers)
}

#' @export
simulate_gym = function(n = 10) {
  exercises = c(rep('Squat', 3), rep('Bicep', 3),
                'Stretching', rep('Bench', 2), rep('Legs', 3),
                'Back', 'Dumbbell', 'Dead lift',
                'Pullup', 'Weighted pullup', 'Hammer')
  days = c('Monday', rep('Tuesday', 2), 'Wednesday',
           'Thursday', 'Rest day', rep('No leg day', 4))
  nums = sample((round(n/2)):(n*3), length(exercises), replace=T)

  gym = data.frame(unname(do.call(c, mapply(function(d, n) {rep(d, n)}, exercises, nums))))
  colnames(gym) = 'exercises'
  gym['days'] = sample(days, nrow(gym), replace=T)
  return(gym)
}

max_stars_repo_name: teng-gao/istack

max_issues_repo_name: teng-gao/istack

max_forks_repo_name: teng-gao/istack


text:
args<-commandArgs(T)

library(RCytoscape)
g <- new('graphNEL', edgemode='undirected')

## read node file
node = read.table(args[1], sep="\t", comment.char = "`")
for (i in 1:nrow(node)){
	g <- graph::addNode(as.character(node[i,1]), g)
}

## read edge file
edge = read.table(args[2], sep="\t")
for (i in 1:nrow(edge)){
	g <- graph::addEdge(as.character(edge[i,1]), as.character(edge[i,2]), g)
}

## establish connection to Cytoscape
cy <- CytoscapeConnection()
deleteAllWindows(cy)
## init
setDefaultBackgroundColor (cy, '#ffffff')
cw <- new.CytoscapeWindow (args[3], graph=g)
setDefaultNodeShape (cw, 'round_rect')
setWindowSize (cw, 1612, 1231)
displayGraph (cw)
## modif node
for (i in 1:nrow(node)){
	setNodeColorDirect (cw, as.character(node[i,1]), as.character(node[i,2]))
}
## adjust network layout
layoutNetwork (cw, 'force-directed')
redraw (cw)
## adjust img size
zoom <- getZoom(cw) * 0.9
setZoom(cw, zoom)
## save img
id = paste(args[3], ".png", sep="")
saveImage (cw, id, 'png', 1.0)


text:
######################################################################
###
###  Decision tree -  Titanic
### The purpose of this dataset is to predict which people
### are more likely to survive after the collision with the iceberg.
### The dataset contains 13 variables and 1309 observations.
### The dataset is ordered by the variable X.
### https://www.guru99.com/r-decision-trees.html#4
###
#######################################################################

###  Variable	Definition	Key
### survival	Survival		0 = No, 1 = Yes
### pclass	Ticket class	1 = 1st, 2 = 2nd, 3 = 3rd
### sex		Sex
### Age		Age in years
### sibsp		# of siblings / spouses aboard the Titanic
### parch		# of parents / children aboard the Titanic
### ticket	Ticket number
### fare		Passenger fare
### cabin		Cabin number
### embarked	Port of Embarkation	C = Cherbourg, Q = Queenstown, S = Southampton
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ###

rm(list=ls())

### Libraries
library(tree)
library(dplyr)
library(rpart.plot)
library(rpart)
library(ROCR)
library(randomForest)
library(caret)
library(e1071)
# Questa funzione installa una librerie che servirà per i random forest
#options(repos='http://cran.rstudio.org')
#have.packages <- installed.packages()
#cran.packages <- c('devtools','plotrix','randomForest','tree')
#to.install <- setdiff(cran.packages, have.packages[,1])
#if(length(to.install)>0) install.packages(to.install)

#library(devtools)
#if(!('reprtree' %in% installed.packages())){
#  install_github('araastat/reprtree')
#}
#for(p in c(cran.packages, 'reprtree')) eval(substitute(library(pkg), list(pkg=p)))
library(reprtree)
### ### ### ### ### ###

set.seed(678)
DIR     = "/Users/gianlucamastrantonio/Desktop/statistico/lavori/lezioni/DataSpace/2018/Lez7/"
Dataset = read.csv(paste(DIR,"titanic.csv",sep=""))
head(Dataset)
tail(Dataset)

# Puliamo il dataset
#
DatasetClean <- Dataset %>%
  select(-c(home.dest, cabin, name, X, ticket)) %>%
  #Convert to factor level
  mutate(pclass = factor(pclass, levels = c(1, 2, 3), labels = c('Upper', 'Middle', 'Lower')),
         survived = factor(survived, levels = c(0, 1), labels = c('No', 'Yes'))) %>%
  na.omit()
glimpse(DatasetClean)
summary(DatasetClean)
# dividiamo il campioni in train and test
shuffle_index = sample(1:nrow(DatasetClean))
head(shuffle_index)
DatasetShuffled = DatasetClean[shuffle_index,]

### Creiamo una funzione per fare lo split
create_train_test  =  function(data, size = 0.8) {
	#
  n_row = nrow(data)
  total_row = round(size*n_row)
  train_sample =   1: total_row
  ret = list(Train = data[train_sample,],Test=data[-train_sample,])
  return(ret)
}
App =  create_train_test(DatasetShuffled, 0.8)
DataTrain = App$Train
DataTest = App$Test
dim(DataTrain)
dim(DataTest)
rm(App)

# controlliamo che la randomizzazione abbia funzionato
prop.table(table(DataTrain$survived))
prop.table(table(DataTest$survived))

# automatizziamo il processo

# versione 1
PropTable = function(x)
{
  prop.table(table(x))
}
apply(DataTrain[,c(1,2,3,5,6,8)],2,PropTable)
apply(DataTest[,c(1,2,3,5,6,8)],2,PropTable)

# versione 2
IndexVar = c(1,2,3,5,6,8)
for(i in IndexVar)
{
  print(colnames(DataTest)[i])
  print(prop.table(table(DataTrain[,i])))
  print(prop.table(table(DataTest[,i])))
}

# il modello
fit = rpart(survived~., data = DataTrain, method = 'class')
rpart.plot(fit, extra = 1,type=4)
rpart.plot(fit, extra = 106,type=4)
rpart.plot(fit, extra = 100,type=4, under=T)
#?rpart.plot  per vedere i vari tipi grafici
rpart.rules(fit)

# facciamo previsione (classe)
Pred = predict(fit, DataTest, type = 'class')
plot(Pred)
# Controlliamo la previsione
Tab =  table(DataTest$survived, Pred)
Tab
Accuracy = sum(diag(Tab))/sum(Tab)

# curva roc
# stimiamo le probs
PredProb = predict(fit, DataTest)
summary(PredProb)
pred = prediction(PredProb[,2], DataTest$survived)
perf = performance(pred,"tpr","fpr")
plot(perf,colorize=TRUE)
abline(a=0,b=1)

## vadiamo meglio come funziona la curca ROC
i = 5
a = perf@alpha.values[[1]][i]
fp = perf@x.values[[1]][i]
tp = perf@y.values[[1]][i]
glm.pred=rep("Down",252)
X = ifelse(PredProb[,2]>=a,"Yes","No")
table(X,DataTest$survived)
57/(37+51)
13/(108+13)

# tuning degli hyperparametri

# creaiamo uan funzione che mostre l'accuratezza
accuracy_tune = function(fit) {
  predict_unseen = predict(fit, DataTest, type = 'class')
  table_mat = table(DataTest$survived, predict_unseen)
  accuracy_Test = sum(diag(table_mat)) / sum(table_mat)
  accuracy_Test
}
# per fare il tuning dobbiamo settare
# rpart.control(minsplit = 20, minbucket = round(20/3), maxdepth = 30)
# tra i molti parametri abbiamo
# -minsplit: Set the minimum number of observations in the node before the algorithm perform a split
# -minbucket:  Set the minimum number of observations in the final note i.e. the leaf
# -maxdepth: Set the maximum depth of any node of the final tree. The root node is treated a depth 0

# facciamo untest e vediamo come migliorare l'accuratezza
control = rpart.control(minsplit = 4,
                         minbucket = round(5 / 3),
                         maxdepth = 3,
                         cp = 0)
tune_fit <- rpart(survived~., data = DataTrain, method = 'class', control = control)
accuracy_tune(tune_fit)
# e vediamo la nuova roc
PredProbTune = predict(tune_fit, DataTest)
summary(PredProbTune)
predTune = prediction(PredProbTune[,2], DataTest$survived)
perfTune = performance(predTune,"tpr","fpr")
plot(perfTune,colorize=TRUE)
abline(a=0,b=1)
plot(perf, add=T,colorize=TRUE, lty=2)

### ### ### ### ### ### ### ### ### ### ###
### Random forest
### ### ### ### ### ### ### ### ### ### ###

# modello base
RF = randomForest(survived~. , ntree=10, mtry=1,data = DataTrain)
plot(RF)
RF
varImpPlot(RF)
# alcuni parametri
# ntree: number of trees in the forest
# mtry: Number of candidates draw to feed the algorithm. By default, it is the square of the number of columns.
# maxnodes: Set the maximum amount of terminal nodes in the forest
# importance=TRUE: Whether independent variables importance in the random forest be assessed

## plottiamo l'albero 2
reprtree:::plot.getTree(RF,k=2)

# Cerchiamo il miglior set di parametri

# troviamo mtry
tuneGrid <- expand.grid(mtry = c(1: 10))

trControl <- trainControl(method = "cv",
                          number = 10,
                          search = "grid")
rf_mtry <- train(survived~.,
                 data = DataTrain,
                 method = "rf",
                 metric = "Accuracy",
                 tuneGrid = tuneGrid,
                 trControl = trControl,
                 importance = TRUE,
                 nodesize = 14,
                 ntree = 300)
rf_mtry

best_mtry <- rf_mtry$bestTune$mtry

# cerchiamo il miglior maxnode
store_maxnode <- list()
tuneGrid <- expand.grid(.mtry = best_mtry)
for (maxnodes in seq(6,30,by=3)) {
  #set.seed(1234)
  rf_maxnode <- train(survived~.,
                      data = DataTrain,
                      method = "rf",
                      metric = "Accuracy",
                      tuneGrid = tuneGrid,
                      trControl = trControl,
                      importance = TRUE,
                      nodesize = 14,
                      maxnodes = maxnodes,
                      ntree = 300)
  current_iteration <- toString(maxnodes)
  store_maxnode[[current_iteration]] <- rf_maxnode
}
results_mtry <- resamples(store_maxnode)
summary(results_mtry)

# cerchiamo il miglior ntree
store_maxtrees <- list()
for (ntree in c(250, 300, 350, 400, 450, 500, 550, 600, 800, 1000, 2000)) {
  set.seed(5678)
  rf_maxtrees <- train(survived~.,
                       data = DataTrain,
                       method = "rf",
                       metric = "Accuracy",
                       tuneGrid = tuneGrid,
                       trControl = trControl,
                       importance = TRUE,
                       nodesize = 14,
                       maxnodes = 24,
                       ntree = ntree)
  key <- toString(ntree)
  store_maxtrees[[key]] <- rf_maxtrees
}
results_tree <- resamples(store_maxtrees)
summary(results_tree)

# best model
fit_rf <- train(survived~.,
                DataTrain,
                method = "rf",
                metric = "Accuracy",
                tuneGrid = tuneGrid,
                trControl = trControl,
                importance = TRUE,
                nodesize = 14,
                ntree = 800,
                maxnodes = 24)

# stimiamo il modello con random Forest
RFtot = randomForest(survived~. , ntree=800, maxnodes = 24,nodesize = 14,mtry=best_mtry,data = DataTrain)
varImp(RFtot)
varImpPlot(RFtot)

## plottiamo l'albero 2
reprtree:::plot.getTree(RFtot,k=2)

# facciamo previsione
prediction <-predict(RFtot, DataTest)

# possiamo vedere la matrice di confusione
confusionMatrix(prediction, DataTest$survived)

# curva ROC
PredProbRF = predict(fit, DataTest,  type="prob")
summary(PredProbRF)
predRF = prediction(PredProbRF[,2], DataTest$survived)
perfRF = performance(pred,"tpr","fpr")
plot(perfRF,colorize=TRUE)
abline(a=0,b=1)
plot(perfTune,colorize=TRUE,add=T)


text:
# Q3
# search atomic

# Q4
x<-4
class(x)

# Q5
y <- c(4, TRUE)
y
class(y)

# Q6
a <- c(3, 2, 10)
b <- c(1,3, 5)
cbind(a,b)

#Q7
dim(a)

# Q8
c <- list(2, "a", "b", TRUE)
c[[1]]

# Q9
d <- 1:4
e <- 2:3
class(d+e)

# Q 10
f <- c(17, 14, 4, 5, 13, 12, 10)
f[f > 10] <- 4

# Q11, 12, 13, 14
names(hw1_data)
nrow(hw1_data)
hw1_data[1:2,]
tail(hw1_data)

## Q 15-20

# What is the value of Ozone in the 47th row?
hw1_data[47,]

# How many missing values are in the Ozone column of this data frame
sum(is.na(hw1_data$Ozone))
# What is the mean of the Ozone column in this dataset? Exclude missing values (coded as NA) from this calculation.
mean(na.omit(hw1_data$Ozone))

# Extract the subset of rows of the data frame where Ozone values are above 31 and Temp values are above 90.
# What is the mean of Solar.R in this subset?
hw1_data.sub1 <- subset(hw1_data, Ozone > 31 & Temp > 90)
mean(hw1_data.sub1$Solar.R)

# What is the mean of "Temp" when "Month" is equal to 6?
hw1_data.sub2 <- subset(na.omit(hw1_data), Month == 6)
mean(hw1_data.sub2$Temp)

# What was the maximum ozone value in the month of May (i.e. Month = 5)?
hw1_data.sub3 <- subset(na.omit(hw1_data), Month == 5)
max(hw1_data.sub3$Ozone)

max_stars_repo_name: rahlk/Rosie

max_issues_repo_name: rahlk/Rosie

max_forks_repo_name: rahlk/Rosie