text:
library(rstan)
library(data.table)
library(lubridate)
library(gdata)
library(dplyr)
library(tidyr)
library(EnvStats)

source('nature/utils/format-data.r')

make_table <- function(filename){

  load(paste0("nature/results/", filename,'-stanfit.Rdata'))
  out <- rstan::extract(fit)
  inf_mean_delta <- vector("list", length = 11)
  inf_lower_delta <- vector("list", length = 11)
  inf_upper_delta <- vector("list", length = 11)

  inf_mean_model <- vector("list", length = 11)
  inf_lower_model <- vector("list", length = 11)
  inf_upper_model <- vector("list", length = 11)

  inf_mean_null <- vector("list", length = 11)
  inf_lower_null <- vector("list", length = 11)
  inf_upper_null <- vector("list", length = 11)

  mean_delta <- vector("list", length = 11)
  lower_delta <- vector("list", length = 11)
  upper_delta <- vector("list", length = 11)

  mean_null <- vector("list", length = 11)
  lower_null <- vector("list", length = 11)
  upper_null <- vector("list", length = 11)

  mean_model <- vector("list", length = 11)
  lower_model <- vector("list", length = 11)
  upper_model <- vector("list", length = 11)

  cf_deaths_overall = rep(0,2000)
  our_deaths_overall = rep(0,2000)
  cf_infection_overall = rep(0,2000)
  our_infection_overall = rep(0,2000)
  for(i in 1:length(countries)){
    fit0 = make_data(filename, "Null model",counterfactual = TRUE, i, days = 0)
    fit1 = make_data(filename, "Our model",counterfactual = FALSE, i, days= 0)

    cf_deaths_overall = cf_deaths_overall + rowSums(fit0$estimated_deaths)
    our_deaths_overall = our_deaths_overall + rowSums(fit1$estimated_deaths)

    cf_infection_overall = cf_infection_overall + rowSums(fit0$estimated_infections)
    our_infection_overall = our_infection_overall + rowSums(fit1$estimated_infections)

    inf.delta =  fit0$estimated_infections - fit1$estimated_infections
    inf.delta.cum = rowSums(inf.delta)
    inf.probs = quantile(inf.delta.cum,probs=c(.025,.5,.975))
    inf_mean_delta[[i]] <- mean(inf.delta.cum)
    inf_lower_delta[[i]] <- inf.probs[1]
    inf_upper_delta[[i]] <- inf.probs[3]

    inf.cum = rowSums(fit0$estimated_infections)
    inf.probs = quantile(inf.cum,probs=c(.025,.5,.975))
    inf_mean_null[[i]] <- mean(inf.cum)
    inf_lower_null[[i]] <- inf.probs[1]
    inf_upper_null[[i]] <- inf.probs[3]

    inf.cum = rowSums(fit1$estimated_infections)
    inf.probs = quantile(inf.cum,probs=c(.025,.5,.975))
    inf_mean_model[[i]] <- mean(inf.cum)
    inf_lower_model[[i]] <- inf.probs[1]
    inf_upper_model[[i]] <- inf.probs[3]

    death.delta =  fit0$estimated_deaths - fit1$estimated_deaths
    death.delta.cum = rowSums(death.delta)
    death.probs = quantile(death.delta.cum,probs=c(.025,.5,.975))
    mean_delta[[i]] <- mean(death.delta.cum)
    lower_delta[[i]] <- death.probs[1]
    upper_delta[[i]] <- death.probs[3]

    death.cum = rowSums(fit0$estimated_deaths)
    death.probs = quantile(death.cum,probs=c(.025,.5,.975))
    mean_null[[i]] <- mean(death.cum)
    lower_null[[i]] <- death.probs[1]
    upper_null[[i]] <- death.probs[3]

    death.cum = rowSums(fit1$estimated_deaths)
    death.probs = quantile(death.cum,probs=c(.025,.5,.975))
    mean_model[[i]] <- mean(death.cum)
    lower_model[[i]] <- death.probs[1]
    upper_model[[i]] <- death.probs[3]
  }
  countries <- c(
    "Denmark",
    "Italy",
    "Germany",
    "Spain",
    "United_Kingdom",
    "France",
    "Norway",
    "Belgium",
    "Austria",
    "Sweden",
    "Switzerland",
    "Total"
  )
  # data frame deaths averted
  df_delta <- data.frame( unlist(mean_delta), unlist(lower_delta),  unlist(upper_delta))
  names(df_delta) <-c('mean', 'lower', 'higher')
  all_q = quantile(cf_deaths_overall-our_deaths_overall,prob=c(.025,.975))
  x <- data.frame(mean(cf_deaths_overall-our_deaths_overall),all_q[1],all_q[2])
  names(x)  <-c('mean', 'lower', 'higher')
  df_delta <- rbind(df_delta,x)
  df_delta <- data.frame(apply(signif(df_delta[,1:3],2),2,prettyNum,big.mark=','))
  df_delta['text'] <- paste0(df_delta[,1]," [",df_delta[,2]," - ",df_delta[,3],"]")
  df_delta <- cbind(df_delta, countries)
  write.csv(df_delta[,c('countries', 'text')], file = paste0('nature/results/deaths-averted-',filename, '.csv'))

  # deaths null
  df_null <- data.frame( unlist(mean_null), unlist(lower_null),  unlist(upper_null))
  names(df_null) <-c('mean', 'lower', 'higher')
  all_q = quantile(cf_deaths_overall,prob=c(.025,.975))
  x <- data.frame(mean(cf_deaths_overall),all_q[1],all_q[2])
  names(x)  <-c('mean', 'lower', 'higher')
  df_null <- rbind(df_null,x)
  df_null <- data.frame(apply(signif(df_null[,1:3],2),2,prettyNum,big.mark=','))
  df_null['text'] <- paste0(df_null[,1]," [",df_null[,2]," - ",df_null[,3],"]")
  df_null <- cbind(df_null, countries)
  write.csv(df_null[,c('countries', 'text')], file = paste0('nature/results/deaths-null-',filename, '.csv'))

  # deaths model
  df_model <- data.frame( unlist(mean_model), unlist(lower_model),  unlist(upper_model))
  names(df_model) <-c('mean', 'lower', 'higher')
  all_q = quantile(our_deaths_overall,prob=c(.025,.975))
  x <- data.frame(mean(our_deaths_overall),all_q[1],all_q[2])
  names(x)  <-c('mean', 'lower', 'higher')
  df_model <- rbind(df_model,x)
  df_model <- data.frame(apply(signif(df_model[,1:3],2),2,prettyNum,big.mark=','))
  df_model['text'] <- paste0(df_model[,1]," [",df_model[,2]," - ",df_model[,3],"]")
  df_model <- cbind(df_model, countries)
  write.csv(df_model[,c('countries', 'text')], file = paste0('nature/results/deaths-model-',filename, '.csv'))

  # cases model
  df_model <- data.frame( unlist(inf_mean_model), unlist(inf_lower_model),  unlist(inf_upper_model))
  names(df_model) <-c('mean', 'lower', 'higher')
  all_q = quantile(our_infection_overall,prob=c(.025,.975))
  x <- data.frame(mean(our_infection_overall),all_q[1],all_q[2])
  names(x)  <-c('mean', 'lower', 'higher')
  df_model <- rbind(df_model,x)
  df_model <- data.frame(apply(signif(df_model[,1:3],2),2,prettyNum,big.mark=','))
  df_model['text'] <- paste0(df_model[,1]," [",df_model[,2]," - ",df_model[,3],"]")
  df_model <- cbind(df_model, countries)
  write.csv(df_model[,c('countries', 'text')], file = paste0('nature/results/cases-model-',filename, '.csv'))

  df_pop = readRDS('nature/data/popt-ifr.rds')
  df_pop = df_pop[1:11, c('country', 'popt')]
  df_pop <- df_pop[match(countries[1:11], df_pop$country),]

  # % cases model
  df_model <- data.frame( unlist(inf_mean_model), unlist(inf_lower_model),  unlist(inf_upper_model))
  df_model <- (df_model / df_pop$popt) * 100
  names(df_model) <-c('mean', 'lower', 'higher')
  all_q = quantile(our_infection_overall,prob=c(.025,.975))
  x <- data.frame(mean(our_infection_overall),all_q[1],all_q[2])
  x <- (x/  sum(df_pop$popt)) * 100
  names(x)  <-c('mean', 'lower', 'higher')
  df_model <- rbind(df_model,x)
  df_model <- data.frame(apply(signif(df_model[,1:3],2),2,prettyNum,big.mark=','))
  df_model['text'] <- paste0(df_model[,1]," [",df_model[,2]," - ",df_model[,3],"]")
  df_model <- cbind(df_model, countries)
  write.csv(df_model[,c('countries', 'text')], file = paste0('nature/results/per-cases-model-',filename, '.csv'))
}


text:
# 2. faza: Uvoz podatkov
source("lib/libraries.r", encoding="UTF-8")
moji_podatki <- read_delim("podatki/goli_podatki.csv",
                           ";", escape_double=FALSE, trim_ws=TRUE,
                           locale=locale(encoding="windows-1250",
                                         decimal_mark=",", grouping_mark="."),
                           col_types=cols(DATUM_UVOZA=col_date(format="%d.%m.%Y"),
                                          `Osnovna sredstva 2018`=col_number()),
                           na=c("#VALUE!", "/", "0", "", "-", "ni podatka"), # ali je 0 res manjkajoč podatek?
                           n_max=632)

#moji_podatki <- read_delim("podatki/goli_podatki.csv",
#                           ";", escape_double = FALSE, trim_ws = TRUE, locale=locale(encoding="windows-1250"),
#                           na=c("#VALUE!","/", "0",""))

#OBDELAVA PODATKOV
moji_podatki$NAZIV <- gsub("\\?", "č", moji_podatki$NAZIV,ignore.case = FALSE)
moji_podatki$NASLOV <- gsub("\\?", "č", moji_podatki$NASLOV,ignore.case = FALSE)
moji_podatki$KRAJ <- gsub("\\?", "č", moji_podatki$KRAJ,ignore.case = FALSE)
moji_podatki$DEJAVNOST <- gsub("\\?", "č", moji_podatki$DEJAVNOST,ignore.case = FALSE)
moji_podatki$VELIKOST_SUBJEKTA <- gsub("\\?", "č", moji_podatki$VELIKOST_SUBJEKTA,ignore.case = FALSE)

osnovni_podatki <- moji_podatki[c(1:632),c(4,8,10,11,13,14)]
stevilski_podatki <- moji_podatki[c(1:632),c(4,8,11,22:29,32,34,36,38,40,41)]

colnames(osnovni_podatki)[2:6] <- c("REGIJA","DEJAVNOST","DEJAVNOST ŠTEVILKA","VELIKOST SUBJEKTA","ŠTEVILO ZAPOSLENIH")
colnames(stevilski_podatki)[6:7] <- c("Profit_2018","Profit_2017")
colnames(stevilski_podatki)[14:15] <- c("Profit_2016","Profit_2015")

#LEGENDA REGIJ
legenda_regij <- matrix(c("Ljubljana","Maribor","Celje","Kranj","Nova Gorica", "Koper", "Novo mesto", "Murska Sobota",
                          1,2,3,4,5,6,8,9),ncol=2)
colnames(legenda_regij) <- c('KRAJ', 'REGIJA')

legenda_regij <- tibble::as_tibble(legenda_regij)
legenda_regij$`REGIJA` <- parse_integer(legenda_regij$`REGIJA`)

#LEGENDA DEJAVNOSTI
legenda_dejavnosti <- osnovni_podatki[c("DEJAVNOST", "DEJAVNOST ŠTEVILKA")]
legenda_dejavnosti <- distinct(legenda_dejavnosti)
#legenda_dejavnosti <- legenda_dejavnosti[arrange(legenda_dejavnosti$`DEJAVNOST ŠTEVILKA`)]

prihodki_regija <- stevilski_podatki[c(1:632),c(2,4,5,12,13)]
#NOVA <-prihodki_regija %>%
  #group_by(regija) %>%
  #summarise(prihodek2018 = sum(`Prihodki 2018`, na.rm = TRUE),
            #prihdek2017 = sum(`Prihodki 2017`, na.rm = TRUE))

#noče mi seštevati po več stolpcih hkrati.
#prihodki_regija$`Prihodki_2017` <- as.numeric(prihodki_regija$`Prihodki_2017`)

prihodki_regija <- gather(prihodki_regija, "leto", "Dohodek", 2:5)
prihodki_regija$Dohodek <- gsub(",.*", "", prihodki_regija$Dohodek)
prihodki_regija$Dohodek <- gsub("\\.", "", prihodki_regija$Dohodek)
prihodki_regija$`Dohodek` <- parse_integer(prihodki_regija$`Dohodek`)
prihodki_regija$leto <- gsub("Prihodki", "", prihodki_regija$leto)
prihodki_regija[is.na(prihodki_regija)] <- 0

# graf prihodki glede na regijo je vredu

dobicek_regija <- stevilski_podatki %>% select(2, 6, 7, 14, 15)
#dobicek_regija$`Profit_2018` <- as.numeric(dobicek_regija$`Profit_2018`)
#dobicek_regija$`Profit_2017` <- as.numeric(dobicek_regija$`Profit_2017`)
dobicek_regija <- gather(dobicek_regija, "leto", "Profit", 2:5)
dobicek_regija$Profit <- gsub(",.*", "", dobicek_regija$Profit)
dobicek_regija$Profit <- gsub("\\.", "", dobicek_regija$Profit)
dobicek_regija$`Profit` <- parse_integer(dobicek_regija$`Profit`)
dobicek_regija$leto <- gsub("Profit_", "", dobicek_regija$leto)
dobicek_regija[is.na(dobicek_regija)] <- 0

# graf dobicek glede na regijo je vredu

zapolseni_regija_dejavnost <- osnovni_podatki %>% select(2, 3, 6)

#tabela prihodkov dobičkov in zaposlenih
prihodki_zaposleni <- inner_join(osnovni_podatki, stevilski_podatki, by = "NAZIV") %>% select(3, 9:10, 17:18, 6)
prihodki_zaposleni <- gather(prihodki_zaposleni, "leto", "Dohodek", 2:5)
prihodki_zaposleni$Dohodek <- gsub(",.*", "", prihodki_zaposleni$Dohodek)
prihodki_zaposleni$Dohodek <- gsub("\\.", "", prihodki_zaposleni$Dohodek)
prihodki_zaposleni$`Dohodek` <- parse_integer(prihodki_zaposleni$`Dohodek`)
prihodki_zaposleni$leto <- gsub("Prihodki", "", prihodki_zaposleni$leto)
prihodki_zaposleni[is.na(prihodki_zaposleni)] <- 0

do10 <- prihodki_zaposleni %>% filter(prihodki_zaposleni$`ŠTEVILO ZAPOSLENIH` < 10)

od10do30 <- prihodki_zaposleni %>% filter(prihodki_zaposleni$`ŠTEVILO ZAPOSLENIH` <= 30, prihodki_zaposleni$`ŠTEVILO ZAPOSLENIH` >= 10)

od30do100 <- prihodki_zaposleni %>% filter(prihodki_zaposleni$`ŠTEVILO ZAPOSLENIH` <= 100, prihodki_zaposleni$`ŠTEVILO ZAPOSLENIH` > 30)

nad100 <- prihodki_zaposleni %>% filter(prihodki_zaposleni$`ŠTEVILO ZAPOSLENIH` > 100)

dobicek_zaposleni <- inner_join(osnovni_podatki, stevilski_podatki, by = "NAZIV") %>% select(3, 11:12, 19:20, 6)
dobicek_zaposleni <- gather(dobicek_zaposleni, "Leto", "Profit", 2:5)
dobicek_zaposleni$Profit <- gsub(",.*", "", dobicek_zaposleni$Profit)
dobicek_zaposleni$Profit <- gsub("\\.", "", dobicek_zaposleni$Profit)
dobicek_zaposleni$`Profit` <- parse_integer(dobicek_zaposleni$`Profit`)
dobicek_zaposleni$Leto <- gsub("Profit_", "", dobicek_zaposleni$Leto)
dobicek_zaposleni[is.na(dobicek_zaposleni)] <- 0
dobicek_zaposleni$DEJAVNOST <- gsub("Zobozdravstvena dejavnost", "Zobozdr. dej.", dobicek_zaposleni$DEJAVNOST)
dobicek_zaposleni$DEJAVNOST <- gsub("Specialistična zunajbolnišnična zdravstvena dejavnost", "Spec. zunajbol. zdr. dej.", dobicek_zaposleni$DEJAVNOST)
dobicek_zaposleni$DEJAVNOST <- gsub("Druge zdravstvene dejavnosti", "Dr. zdr. dej.", dobicek_zaposleni$DEJAVNOST)
dobicek_zaposleni$DEJAVNOST <- gsub("Dejavnost nastanitvenih ustanov za oskrbo starejših in invalidnih oseb", "Dej. nast. ustan. za star.", dobicek_zaposleni$DEJAVNOST)
dobicek_zaposleni$DEJAVNOST <- gsub("Bolnišnična zdravstvena dejavnost", "Bol. zdr. dej.", dobicek_zaposleni$DEJAVNOST)
dobicek_zaposleni$DEJAVNOST <- gsub("Splošna zunajbolnišnična zdravstvena dejavnost", "Spl. zunajbol. zdr.", dobicek_zaposleni$DEJAVNOST)
dobicek_zaposleni$DEJAVNOST <- gsub("Dejavnost nastanitvenih ustanov za oskrbo duševno prizadetih, duševno obolelih in zasvojenih oseb", "Dej. nast. ustan. za duš. priz.", dobicek_zaposleni$DEJAVNOST)
dobicek_zaposleni$DEJAVNOST <- gsub("Alternativne oblike zdravljenja", "Alt. obl. zdr.", dobicek_zaposleni$DEJAVNOST)
dobicek_zaposleni$DEJAVNOST <- gsub("Dejavnost nastanitvenih ustanov za bolniško nego", "Dej. nast. ustan. za nego", dobicek_zaposleni$DEJAVNOST)

Profit_do10 <- dobicek_zaposleni %>% filter(dobicek_zaposleni$`ŠTEVILO ZAPOSLENIH` < 10)

Profit_od10do30 <- dobicek_zaposleni %>% filter(dobicek_zaposleni$`ŠTEVILO ZAPOSLENIH` <= 30, dobicek_zaposleni$`ŠTEVILO ZAPOSLENIH` >= 10)

Profit_od30do100 <- dobicek_zaposleni %>% filter(dobicek_zaposleni$`ŠTEVILO ZAPOSLENIH` <= 100, dobicek_zaposleni$`ŠTEVILO ZAPOSLENIH` > 30)

Profit_nad100 <- dobicek_zaposleni %>% filter(dobicek_zaposleni$`ŠTEVILO ZAPOSLENIH` > 100)

#

ext: r


text: sum(scan("", numeric(0), 2))

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
# Load packages
require(ggplot2)
require(grid)
require(gridExtra)
require(RColorBrewer)
require(dplyr)
require(tidyr)
require(stringr)
require(xtable)
require(GGally)
require(portfolio)
require(kableExtra)
require(captioner)
tables <- captioner(prefix = "Tabell")
figures <- captioner(prefix = "Figur")

# Sourcing the hietala_theme
source("hietala_theme.R")

# Loading data used for examples
data(iris)

exempeldata <- read.csv2(file = "data_sets/732G45_exempeldata.csv")

max_stars_repo_name: hietalai/hietala-website

max_issues_repo_name: hietalai/hietala-website

max_forks_repo_name: hietalai/hietala-website


text:
#' @title whichEvents
#' @description unknown
#' @param \code{lonlatpoint} is the central data point you want to find the surroundings c(lon,lat)
#' @param \code{radius} is the radius of the buffer zone you want to identify
#' @param \code{event.data} is the data frame for finding the points within a specified radius
#' @param \code{on.boundary} if true will include all points on the boundary of the polygon
#' @family abysmally documented
#' @author  unknown, \email{<unknown>@@dfo-mpo.gc.ca}
#' @export
whichEvents <- function(lonlatpoint,radius,event.data,on.boundary=T) {
	loadPackages(PBSmapping)

	poly <- makePBS(bufferCircle(lonlat=lonlatpoint,radius=radius),polygon=T)
	e <- makePBS(event.data,polygon=F)
	e <- e[e$EID %in% findPolys(e,poly)[,1],]

	oo <- list(lonlatpoint,radius,poly,EventsinPolygon=e)

	return(oo)
}


text:
> permutation(letters[1:3])
[[1]]
[1] "c" "b" "a"

[[2]]
[1] "b" "c" "a"

[[3]]
[1] "b" "a" "c"

[[4]]
[1] "c" "a" "b"

[[5]]
[1] "a" "c" "b"

[[6]]
[1] "a" "b" "c"

max_stars_repo_name: mullikine/RosettaCodeData


text:
# Goal: Show the efficiency of the mean when compared with the median
#       using a large simulation where both estimators are applied on
#       a sample of U(0,1) uniformly distributed random numbers.

one.simulation <- function(N=100) {     # N defaults to 100 if not supplied
  x <- runif(N)
  return(c(mean(x), median(x)))
}

# Simulation --
results <- replicate(100000, one.simulation(20)) # Gives back a 2x100000 matrix

# Two kernel densities --
k1 <- density(results[1,])              # results[1,] is the 1st row
k2 <- density(results[2,])

# A pretty picture --
xrange <- range(k1$x, k2$x)
plot(k1$x, k1$y, xlim=xrange, type="l", xlab="Estimated value", ylab="")
grid()
lines(k2$x, k2$y, col="red")
abline(v=.5)
legend(x="topleft", bty="n",
       lty=c(1,1),
       col=c("black", "red"),
       legend=c("Mean", "Median"))

max_stars_repo_path: RFrontEndSolution/Tests Repo/Rscripts All (163 files)/p1.r


text:
library(tidyverse)
library(readr)
library(ggthemes)

xyj_result <- read_csv("python_study/xyj/result.csv")
View(xyj_result)

Top10 <- xyj_result[1:10,]
p = ggplot(Top10, aes(x = Word, y = Number, fill = Word))+geom_bar(stat = "identity") + theme(text = element_text(family = 'SimSun'))
p


text:
# 0_initialize.r
#
# Created 20190319

library(arm)
library(broom)
library(broom.mixed)
library(car)
library(conflicted)
library(corrplot)
library(feather)
library(geepack)
library(ggExtra)
library(ggstatsplot)
library(ggplot2)
library(gridExtra)
library(here)
library(lme4)
library(lubridate)
library(paletteer)
library(performance)
library(readxl)
library(skimr)
library(tidyr)
library(tidyverse)
library(varhandle)

set.seed(1)
conflict_prefer('here', 'here')

source(here('r', '0_set_up_environment.r'))
source(here('r', '1_import.r'))
source(here('r', '2_preprocess.r'))
source(here('r', '3_model.r'))

#source(here('r', '3_fit_models.r'))
#source(here('r', '4_export.r'))

message('Run successfully.')


text:
There are 17 section headers, starting at offset .*:

Section Headers:
  \[Nr\] Name              Type            Addr     Off    Size   ES Flg Lk Inf Al
  \[ 0\]                   NULL            00000000 000000 000000 00      0   0  0
  \[ 1\] \.hash             HASH            00008000 001000 000038 04   A  2   0  4
  \[ 2\] \.dynsym           DYNSYM          00008038 001038 000090 10   A  3   6  4
  \[ 3\] \.dynstr           STRTAB          000080c8 0010c8 00001d 00   A  0   0  1
  \[ 4\] \.rela\.got         RELA            000080e8 0010e8 000018 0c  AI  2  10  4
  \[ 5\] \.rela\.bss         RELA            00008100 001100 00000c 0c  AI  2  12  4
  \[ 6\] \.dynamic          DYNAMIC         0000810c 00110c 0000b8 08  WA  3   0  4
  \[ 7\] \.rela\.plt         RELA            10000000 002000 00000c 0c  AI  2  10  4
  \[ 8\] \.plt              PROGBITS        10000020 002020 000040 18  AX  0   0 32
  \[ 9\] \.text             PROGBITS        10000060 002060 000040 00  AX  0   0 32
  \[10\] \.got              PROGBITS        100000a0 0020a0 000020 00  WA  0   0  4
  \[11\] \.neardata         PROGBITS        100000c0 0020c0 000004 00  WA  0   0  4
  \[12\] \.bss              NOBITS          100000c4 0020c4 000004 00  WA  0   0  4
  \[13\] \.c6xabi\.attributes C6000_ATTRIBUTES 00000000 0020c4 000019 00      0   0  1
  \[14\] \.symtab           SYMTAB          00000000 0020e0 000170 10     15  20  4
  \[15\] \.strtab           STRTAB          00000000 002250 000057 00      0   0  1
  \[16\] \.shstrtab         STRTAB          00000000 0022a7 00007b 00      0   0  1
Key to Flags:
#...

Elf file type is EXEC \(Executable file\)
Entry point 0x10000060
There are 4 program headers, starting at offset 52

Program Headers:
  Type           Offset   VirtAddr   PhysAddr   FileSiz MemSiz  Flg Align
  LOAD           0x001000 0x00008000 0x00008000 0x001c4 0x001c4 RW  0x1000
  LOAD           0x002000 0x10000000 0x10000000 0x000c4 0x000c8 RWE 0x1000
  DYNAMIC        0x00110c 0x0000810c 0x0000810c 0x000b8 0x000b8 RW  0x4
  GNU_STACK      0x000000 0x00000000 0x00000000 0x00000 0x20000 RWE 0x8

 Section to Segment mapping:
  Segment Sections\.\.\.
   00     \.hash \.dynsym \.dynstr \.rela\.got \.rela\.bss \.dynamic
   01     \.rela\.plt \.plt \.text \.got \.neardata \.bss
   02     \.dynamic
   03

Dynamic section at offset 0x110c contains 18 entries:
  Tag        Type                         Name/Value
 0x00000001 \(NEEDED\)                     Shared library: \[tmpdir/libtestr\.so\]
 0x00000004 \(HASH\)                       0x8000
 0x00000005 \(STRTAB\)                     0x80c8
 0x00000006 \(SYMTAB\)                     0x8038
 0x0000000a \(STRSZ\)                      29 \(bytes\)
 0x0000000b \(SYMENT\)                     16 \(bytes\)
 0x00000015 \(DEBUG\)                      0x0
 0x00000003 \(PLTGOT\)                     0x100000ac
 0x00000002 \(PLTRELSZ\)                   12 \(bytes\)
 0x00000014 \(PLTREL\)                     RELA
 0x00000017 \(JMPREL\)                     0x10000000
 0x00000007 \(RELA\)                       0x80e8
 0x00000008 \(RELASZ\)                     48 \(bytes\)
 0x00000009 \(RELAENT\)                    12 \(bytes\)
 0x70000000 \(C6000_DSBT_BASE\)            0x100000a0
 0x70000001 \(C6000_DSBT_SIZE\)            0x3
 0x70000003 \(C6000_DSBT_INDEX\)           0x0
 0x00000000 \(NULL\)                       0x0

Relocation section '\.rela\.got' at offset 0x10e8 contains 2 entries:
 Offset     Info    Type                Sym\. Value  Symbol's Name \+ Addend
100000b8  00000601 R_C6000_ABS32          100000c0   b \+ 0
100000bc  00000801 R_C6000_ABS32          100000c4   a \+ 0

Relocation section '\.rela\.bss' at offset 0x1100 contains 1 entry:
 Offset     Info    Type                Sym\. Value  Symbol's Name \+ Addend
100000c4  0000081a R_C6000_COPY           100000c4   a \+ 0

Relocation section '\.rela\.plt' at offset 0x2000 contains 1 entry:
 Offset     Info    Type                Sym\. Value  Symbol's Name \+ Addend
100000b4  0000071b R_C6000_JUMP_SLOT      00000000   sub0 \+ 0

Symbol table '\.dynsym' contains 9 entries:
   Num:    Value  Size Type    Bind   Vis      Ndx Name
     0: 00000000     0 NOTYPE  LOCAL  DEFAULT  UND
     1: 10000020     0 SECTION LOCAL  DEFAULT    8
     2: 10000060     0 SECTION LOCAL  DEFAULT    9
     3: 100000a0     0 SECTION LOCAL  DEFAULT   10
     4: 100000c0     0 SECTION LOCAL  DEFAULT   11
     5: 100000c4     0 SECTION LOCAL  DEFAULT   12
     6: 100000c0     4 OBJECT  GLOBAL DEFAULT   11 b
     7: 00000000     0 FUNC    GLOBAL DEFAULT  UND sub0
     8: 100000c4     4 OBJECT  GLOBAL DEFAULT   12 a

Symbol table '\.symtab' contains 23 entries:
   Num:    Value  Size Type    Bind   Vis      Ndx Name
     0: 00000000     0 NOTYPE  LOCAL  DEFAULT  UND
     1: 00008000     0 SECTION LOCAL  DEFAULT    1
     2: 00008038     0 SECTION LOCAL  DEFAULT    2
     3: 000080c8     0 SECTION LOCAL  DEFAULT    3
     4: 000080e8     0 SECTION LOCAL  DEFAULT    4
     5: 00008100     0 SECTION LOCAL  DEFAULT    5
     6: 0000810c     0 SECTION LOCAL  DEFAULT    6
     7: 10000000     0 SECTION LOCAL  DEFAULT    7
     8: 10000020     0 SECTION LOCAL  DEFAULT    8
     9: 10000060     0 SECTION LOCAL  DEFAULT    9
    10: 100000a0     0 SECTION LOCAL  DEFAULT   10
    11: 100000c0     0 SECTION LOCAL  DEFAULT   11
    12: 100000c4     0 SECTION LOCAL  DEFAULT   12
    13: 00000000     0 SECTION LOCAL  DEFAULT   13
    14: 00000000     0 FILE    LOCAL  DEFAULT  ABS tmpdir/shlib-app-1r\.o
    15: 10000060     0 NOTYPE  LOCAL  DEFAULT    9 fish
    16: 00000000     0 FILE    LOCAL  DEFAULT  ABS
    17: 0000810c     0 OBJECT  LOCAL  DEFAULT    6 _DYNAMIC
    18: 100000ac     0 OBJECT  LOCAL  DEFAULT   10 _GLOBAL_OFFSET_TABLE_
    19: 100000a0     0 NOTYPE  LOCAL  DEFAULT   10 __c6xabi_DSBT_BASE
    20: 100000c0     4 OBJECT  GLOBAL DEFAULT   11 b
    21: 00000000     0 FUNC    GLOBAL DEFAULT  UND sub0
    22: 100000c4     4 OBJECT  GLOBAL DEFAULT   12 a

ext: rd


text:
REBOL [
  Title:   "Builds a set of Red Tests to run on an ARM host"
	File: 	 %build-arm-tests.r
	Author:  "Peter W A Wood"
	Version: 0.1.0
	License: "BSD-3 - https://github.com/dockimbel/Red/blob/master/BSD-3-License.txt"
]

;; This script must be run from the Red/red/tests dir

;; supress script messages
store-quiet-mode: system/options/quiet
system/options/quiet: true

;; init
file-chars: charset [#"a" - #"z" #"A" - #"Z" #"0" - #"9" "-" "/"]
a-file-name: ["%" some file-chars ".red" ]
a-test-file: ["--run-test-file-quiet " copy file a-file-name]

target: ask {
Choose ARM target:
1) Linux
2) Android
=> }
target: pick ["Linux-ARM" "Android"] to-integer target

;; make the Arm dir if needed
arm-dir: %runnable/arm-tests/
make-dir/deep arm-dir

;; empty the Arm dir
foreach file read arm-dir [delete join arm-dir file]

;; get the list of test source files
test-files: copy []
all-tests: read %run-all.r
parse/all all-tests [
	thru "Red Units tests"
	any [a-test-file (append test-files to file! file) | skip]
]

;; compile the tests into to runnable/arm-tests

foreach test-file test-files [
  do/args %../../red.r rejoin ["-t " target " " test-file]
  exe: copy find/last/tail test-file "/"
  exe: replace exe ".red" ""
  write/binary join %runnable/arm-tests/ exe read/binary exe
  delete exe
]

;; copy the bash script and mark it as executable
runner: %runnable/arm-tests/run-all.sh
write/binary runner read/binary %run-all.sh
if system/version/4 <> 3 [
	set-modes runner [
	  owner-execute: true
	  group-execute: true
	  world-execute: true
	]
]

;; tidy up
system/options/quiet: store-quiet-mode

print "ARM tests built"


text:
-- @@stderr --
dtrace: failed to compile script test/unittest/typedef/err.D_SYNTAX.TypedefInClause.d: [D_SYNTAX] line 20: syntax error near "typedef"

max_stars_repo_name: alan-maguire/dtrace-utils

max_issues_repo_name: tjfontaine/dtrace-utils

max_forks_repo_name: tjfontaine/dtrace-utils


text:
/*	XOPTypes.r -- Type declaration for XOP specific resources

	This file is used for Macintosh XOPs, not for Windows XOPs.
*/

#include "XOPResources.h"		// Load symbols.

/* XOPI Template -- various information about XOP. */
type 'XOPI' {
	integer;							// XOP protocol version. Use XOP_VERSION.
	integer;							// Code for development system used to make XOP. Use DEV_SYS_CODE.
	integer;							// Obsolete. Set to zero.
	integer;							// Obsolete. Set to zero.
	integer;							// XOP Toolkit version. 600 means version 6.00. Use XOP_TOOLKIT_VERSION. Old XOPs will have zero in this field.
};

/* XOPC Template -- information about XOP command line operations. */
type 'XOPC' {
	integer = $$Countof(CmdArray);
	array CmdArray {
		pstring[31];					/* Operation name -- 31 chars max. */
		integer;						/* Operation's category (see #defines below). */
	};
};

/* XMN1 Template -- description of menus added by XOP. */
type 'XMN1' {
	array XMN1Array {
		integer;						/* Resource ID of 'MENU' resource in XOP. */
		integer;						/* Menu flags. See XOPResources.h. */
	};
};

/* XSM1 Template -- description of submenus added by XOP. */
type 'XSM1' {
	array XSM1Array {
		integer;						/* Resource ID of 'MENU' resource in XOP for submenu. */
		integer;						/* Resource ID of 'MENU' resource in XOP for main menu. */
		integer;						/* Number of item in main menu where submenu is to be attached. */
	};
};

/* XMI1 Template -- description of menu items added by XOP to built-in Igor menus. */
type 'XMI1' {
	array XMI1Array {
		integer;						/* Menu ID of built-in Igor menu. See IgorXOP.h for menu IDs. */
		pstring[31];					/* Text for item to be added to menu. */
		integer;						/* Resource ID of 'MENU' resource in XOP for submenu. */
										/* To be attached to menu item or 0 for no submenu. */
		integer;						/* Menu item flags. See XOPResources.h. */
	};
};

/* XOPF Template -- description of functions added by XOP */
type 'XOPF' {
	integer = $$Countof(XFuncArray);	/* Number of functions added by this XOP. */
	array XFuncArray {
		pstring[31];					/* Function name. */
		integer;						/* Function category. See XOPResources.h. */
		integer;						/* Return value type. See XOPResources.h. */
		integer = $$Countof(ParmArray);	/* Number of parameters to this function. */
		array ParmArray {
			integer;					/* Type of this parameter. See XOPResources.h. */
		};
	};
};

max_stars_repo_name: robnagler/Radia

max_issues_repo_name: loveinthebone/Radia

max_forks_repo_name: loveinthebone/Radia


text:
SPFA <- function(paths=NULL,mdir=NULL,de=NULL,all=NULL,normal=NULL,cancer=NULL){
  nB <- 2000
  pb <- NULL
  pl <- NULL
  leafscores <- NULL
  Names <- NULL
  pathwayIDs <- NULL
  for(i in 1:length(paths)){
    mapkpathway <- try(parseKGML(paste(mdir,paths[i],sep = "/")),TRUE)
    nodes <- nodes(mapkpathway)
    mapkG <- KEGGpathway2Graph(mapkpathway,genesOnly=FALSE,expandGenes=FALSE)
    mapkG1 <- KEGGpathway2Graph(mapkpathway,expandGenes=TRUE)
    mapkG2 <- KEGGpathway2Graph(mapkpathway,genesOnly=FALSE,expandGenes=TRUE)
    allgenes <- nodes(mapkG1)
    allnodes <- nodes(mapkG2)
    allgeneslist <- strsplit(allgenes,"hsa:")
    allgeneslist <- do.call(rbind,allgeneslist)
    allgenes <- allgeneslist[,2]
    allnodeslist <- strsplit(allnodes,"hsa:")
    allnodeslist <- do.call(rbind,allnodeslist)
    allnodes <- allnodeslist[,2]
    g <- igraph.from.graphNEL(mapkG)
    outdegree <- igraph::degree(g,mode ="out")
    Ue <- names(outdegree[which(outdegree==0)])
    Uenodes <- nodes[Ue]
    types <- sapply(Uenodes, getType)
    endnodes <- names(types[which(types=="gene")])
    nonodes <- setdiff(Ue,endnodes)
    indegree<-igraph::degree(g,mode ="in")
    Us <- names(indegree[which(indegree==0)])
    nonodes <- setdiff(nonodes,Us)

    dismatrix <- distances(g,mode="in")
    dismatrix[which(dismatrix=="Inf")] <- 0
    nn <- max(dismatrix)
    dismatrix[which(dismatrix!=0)] <- 1

    ni <- 1
    while(length(nonodes)!=0 & ni<=nn){
      adjnodes<-NULL
      for(j in 1:length(nonodes)){
        adjnode <- neighbors(g,nonodes[j],"in")
        for(ij in 1:length(adjnode)){
          adjtype <- nodes[[adjnode[ij]]]
          nonodes <- setdiff(nonodes,names(adjnode[ij]))
          if(adjtype@type == "gene"){
            endnodes <- c(endnodes,names(adjnode[ij]))
          }
          if(adjtype@type != "gene"){
            adjnodes <- c(adjnodes,names(adjnode[ij]))
          }
        }
      }
      adjnodes <- setdiff(adjnodes,Us)
      nonodes <- adjnodes
      ni <- ni+1
    }
    endnodes <- endnodes[!duplicated(endnodes)]

    endtypes <- nodes[endnodes]
    types1 <- sapply(endtypes, getType)
    endnodes <- names(types1[which(types1=="gene")])
    endgenes <- NULL
    for(ii in 1:length(endtypes)){
      if(length(endtypes[[ii]]) != 0){
        endgene <- endtypes[[ii]]@name
        endgenelist <- strsplit(endgene,"hsa:")
        endgenelist <- do.call(rbind,endgenelist)
        endgene <- endgenelist[,2]
        endgenes <- c(endgenes,endgene)
      }
    }
    endgenes <- endgenes[!duplicated(endgenes)]

    nde <- intersect(de,allgenes)
    pl[i] <- phyper(q=length(nde)-1,m=length(allgenes),n=length(all)-length(allgenes),k=length(de),lower.tail=FALSE)
    conames <- allgenes
    difnames <- setdiff(conames,colnames(normal))
    conames <- intersect(conames,colnames(normal))

    nfor <- normal[,conames]
    nfor <- as.data.frame(nfor)
    cfor <- cancer[,conames]
    cfor <- as.data.frame(cfor)
    Inormal <- cor(nfor)
    Icancer <- cor(cfor)
    Inormal[is.na(Inormal)] <- 1
    Icancer[is.na(Icancer)] <- 1
    Iscore <- Icancer-Inormal
    if(length(difnames)!=0){
      ainter <- rep(0,ncol(Iscore))
      ai <- matrix(rep(ainter,length(difnames)),ncol=length(difnames))
      binter <- rep(0,length(allgenes))
      bi <- matrix(rep(binter,length(difnames)),nrow=length(difnames))
      Iscore <- cbind(Iscore,ai)
      Iscore <- rbind(Iscore,bi)
      genenames <- c(conames,difnames)
      colnames(Iscore) <- genenames
      rownames(Iscore) <- genenames
      Iscore <- Iscore[allgenes,allgenes]
    }
    Iscore <- abs(Iscore)
    notends <- setdiff(allnodes,endgenes)
    g1 <- igraph.from.graphNEL(mapkG2)
    dismatrix1 <- distances(g1,mode="in")
    dismatrix1[which(dismatrix1=="Inf")] <- 0
    dism <- dismatrix1
    colnames(dism) <- allnodes
    rownames(dism) <- allnodes
    dism <- dism[allgenes,allgenes]
    dismatrix1[which(dismatrix1!=0)] <- 1
    colnames(dismatrix1) <- allnodes
    rownames(dismatrix1) <- allnodes
    dismatrix1[notends,] <- 0
    allgenes <- as.character(allgenes)
    dismatrix1 <- dismatrix1[allgenes,allgenes]

    edgesL <- Iscore*dismatrix1/dism
    edgesL[is.na(edgesL)] <- 0
    edgescore <- sum(edgesL)

    zgenes <- NULL
    zgenes1 <- NULL
    for(iiii in 1:ncol(Iscore)){
      if(sum(Iscore[iiii,]) == 0){
        zgenes[iiii] <- rownames(Iscore)[iiii]
      }
      if(sum(Iscore[,iiii]) == 0){
        zgenes1[iiii] <-colnames(Iscore)[iiii]
      }
    }
    zgenes <- zgenes[!is.na(zgenes)]
    zgenes1 <- zgenes1[!is.na(zgenes1)]

    score <- NULL
    for(js in 1:nB){
      esams <- sample(1:ncol(normal),length(allgenes))
      nforsam <- normal[,esams]
      nforsam <- as.data.frame(nforsam)
      cforsam <- cancer[,esams]
      cforsam <- as.data.frame(cforsam)
      Inormalsam <- cor(nforsam)
      Icancersam <- cor(cforsam)
      Inormalsam[is.na(Inormalsam)] <- 1
      Icancersam[is.na(Icancersam)] <- 1
      Iscoresam <- Icancersam-Inormalsam
      colnames(Iscoresam) <- allgenes
      rownames(Iscoresam) <- allgenes
      Iscoresam[zgenes,] <- 0
      Iscoresam[,zgenes1] <- 0
      #Iscoresam <- abs(Iscoresam)*dismatrix1*LLm/dism
      Iscoresam <- abs(Iscoresam)*dismatrix1/dism
      Iscoresam[is.na(Iscoresam)] <- 0
      edgesscore <- sum(Iscoresam)
      score[js] <- edgesscore
    }

    if(edgescore > 0){
      pb[i] <- sum(score >= edgescore)/length(score)
      if(pb[i] <= 0){pb[i]<-1/nB/100}
      if(pb[i] > 1){pb[i]<-1}
    }
    if(edgescore == 0){
      if(all(score == 0)){    #there is nothing to learn from perturbations
        pb[i] <- NA
      }else{
        pb[i] <- 1
      }
    }
    if(edgescore < 0){
      pb[i] <- NA
    }
    Names[i] <- mapkpathway@pathwayInfo@title
    pathwayIDs[i] <- mapkpathway@pathwayInfo@number
    cat("\n");
    cat(paste("Done the ",i," th pathway: ",mapkpathway@pathwayInfo@title,sep=""))
  }
  c <- pb*pl
  pc <- c-c*log(c)
  pcFDR = p.adjust(pc,"fdr")
  pcfwer = p.adjust(pc,"bonferroni")
  res <- data.frame(Names,ID=pathwayIDs,pl=pl,pb=pb,p=pc,pFdr=pcFDR,pFWER=pcfwer,stringsAsFactors=FALSE)

  res <- res[!is.na(res$p),]
  res <- res[order(res$p),]
  rownames(res) <- NULL
  return(res)
}

SPFA_per <- function(mdir = NULL, pathfilename = NULL, normal = NULL, cancer = NULL ){
  mapkpathway<-try(parseKGML(paste(mdir,pathfilename,sep = "/")),TRUE)
  nodes<-nodes(mapkpathway)
  mapkG <- KEGGpathway2Graph(mapkpathway,genesOnly=FALSE,expandGenes=FALSE)
  mapkG1 <- KEGGpathway2Graph(mapkpathway,expandGenes=TRUE)
  mapkG2 <- KEGGpathway2Graph(mapkpathway,genesOnly=FALSE,expandGenes=TRUE)
  allgenes <- nodes(mapkG1)
  allnodes <- nodes(mapkG2)
  allgeneslist <- strsplit(allgenes,"hsa:")
  allgeneslist <- do.call(rbind,allgeneslist)
  allgenes <- allgeneslist[,2]
  allnodeslist <- strsplit(allnodes,"hsa:")
  allnodeslist <- do.call(rbind,allnodeslist)
  allnodes <- allnodeslist[,2]
  g<-igraph.from.graphNEL(mapkG)
  outdegree<-igraph::degree(g,mode ="out")
  Ue<-names(outdegree[which(outdegree==0)])
  Uenodes <- nodes[Ue]
  types <- sapply(Uenodes, getType)
  endnodes <- names(types[which(types=="gene")])
  nonodes <- setdiff(Ue,endnodes)
  indegree<-igraph::degree(g,mode ="in")
  Us<-names(indegree[which(indegree==0)])
  nonodes <- setdiff(nonodes,Us)

  g2 <- igraph.from.graphNEL(mapkG1)
  outdegree1 <- igraph::degree(g2,mode ="out")
  Ue1 <- names(outdegree1[which(outdegree1==0)])
  indegree1 <- igraph::degree(g2,mode ="in")
  Us1 <- names(indegree1[which(indegree1==0)])
  onenodes <- intersect(Ue1,Us1)
  onenodeslist <- strsplit(onenodes,"hsa:")
  onenodeslist <- do.call(rbind,onenodeslist)
  onenodes <- onenodeslist[,2]

  dismatrix<-distances(g,mode="in")
  dismatrix[which(dismatrix=="Inf")]<-0
  nn <- max(dismatrix)
  dismatrix[which(dismatrix!=0)] <- 1

  ni <- 1
  while(length(nonodes)!=0 & ni<=nn){
    adjnodes<-NULL
    for(j in 1:length(nonodes)){
      adjnode <- neighbors(g,nonodes[j],"in")
      for(ij in 1:length(adjnode)){
        adjtype <- nodes[[adjnode[ij]]]
        nonodes <- setdiff(nonodes,names(adjnode[ij]))
        if(adjtype@type=="gene"){
          endnodes <- c(endnodes,names(adjnode[ij]))
        }
        if(adjtype@type!="gene"){
          adjnodes <- c(adjnodes,names(adjnode[ij]))
        }
      }
    }
    adjnodes <- setdiff(adjnodes,Us)
    nonodes <- adjnodes
    ni <- ni+1
  }
  endnodes <- endnodes[!duplicated(endnodes)]
  conames <- allgenes
  difnames <- setdiff(conames,colnames(normal))
  conames <- intersect(conames,colnames(normal))

  endtypes <- nodes[endnodes]
  types1 <- sapply(endtypes, getType)
  endnodes <- names(types1[which(types1=="gene")])
  endgenes <- NULL
  for(ii in 1:length(endtypes)){
    if(length(endtypes[[ii]])!=0){
      endgene <- endtypes[[ii]]@name
      endgenelist <- strsplit(endgene,"hsa:")
      endgenelist <- do.call(rbind,endgenelist)
      endgene <- endgenelist[,2]
      endgenes <- c(endgenes,endgene)
    }
  }

  endgenes <- endgenes[!duplicated(endgenes)]

  nfor <- normal[,conames]
  nfor <- as.data.frame(nfor)
  cfor <- cancer[,conames]
  cfor <- as.data.frame(cfor)
  Inormal <- cor(nfor)
  Icancer <- cor(cfor)
  Inormal[is.na(Inormal)] <- 1
  Icancer[is.na(Icancer)] <- 1
  Iscore <- Icancer-Inormal
  if(length(difnames)!=0){
    ainter <- rep(0,ncol(Iscore))
    ai <- matrix(rep(ainter,length(difnames)),ncol=length(difnames))
    binter <- rep(0,length(allgenes))
    bi <- matrix(rep(binter,length(difnames)),nrow=length(difnames))
    Iscore <- cbind(Iscore,ai)
    Iscore <- rbind(Iscore,bi)
    genenames <- c(conames,difnames)
    colnames(Iscore) <- genenames
    rownames(Iscore) <- genenames
    Iscore <- Iscore[allgenes,allgenes]
  }
  Iscore <- abs(Iscore)
  notends <- setdiff(allnodes,endgenes)
  g1<-igraph.from.graphNEL(mapkG2)
  dismatrix1<-distances(g1,mode="in")
  dismatrix1[which(dismatrix1=="Inf")]<-0
  dism <- dismatrix1
  colnames(dism) <- allnodes
  rownames(dism) <- allnodes
  dism <- dism[allgenes,allgenes]
  dismatrix1[which(dismatrix1!=0)] <- 1
  colnames(dismatrix1) <- allnodes
  rownames(dismatrix1) <- allnodes
  dismatrix1[notends,] <- 0
  allgenes <- as.character(allgenes)
  dismatrix1 <- dismatrix1[allgenes,allgenes]

  edgesL <- Iscore*dismatrix1/dism
  edgesL[is.na(edgesL)] <- 0

  endscore <- apply(edgesL,1,sum)
  endscore <- endscore[which(endscore!=0)]
  endgenes <- setdiff(names(endscore),onenodes)
  endscore <- endscore[endgenes]
  endscore <- -sort(-endscore)

  genes <- factor(names(endscore),levels=names(endscore))
  scores <- data.frame(genes,perturbation=as.numeric(endscore))
  p<-ggplot(scores,aes(x=scores$genes,y=scores$perturbation))+geom_bar(stat="identity")+theme(axis.text.x = element_text(angle = 45, hjust = 0.5, vjust = 0.5))+theme(panel.grid=element_blank())+theme(panel.grid.major = element_line(colour = NA))+theme(legend.position="none")
  return(p+xlab("Effector genes")+ylab("Signal variation received by effector genes"))
}

max_stars_repo_path: SPFA Function.r

max_stars_repo_name: ZhenshenBao/SPFA

max_issues_repo_path: SPFA Function.r

max_issues_repo_name: ZhenshenBao/SPFA

max_forks_repo_path: SPFA Function.r

max_forks_repo_name: ZhenshenBao/SPFA


text:
##
## Fit GEE models
##
library(geepack)
library(MuMIn)
library(arm)

source("data functions/standardize-gee.r")
source("data functions/confint-gee.r")
source("data functions/vcov.averaging.geeglm.r")

source("init.r")

hh.all$Ego.VillageID = factor(hh.all$Ego.VillageID)

##
## intercept-only model
##
gee.null = geeglm(HelpObserved ~ 1, id=Ego.HH, family="binomial", data=hh.all, corstr="exchangeable", scale.fix=T)

##
###
##
gee.village = geeglm(HelpObserved ~ Ego.VillageID, id=Ego.HH, family="binomial", data=hh.all, corstr="exchangeable", scale.fix=T)

##
## control model (distance + size + wealth)
##
gee.control = geeglm(HelpObserved ~ Distance + RelativeWealthRank + Ego.Size + Alter.Size + Ego.VillageID,
                     id=Ego.HH, family="binomial", data=hh.all, corstr="exchangeable", scale.fix=T)

##
## reciprocated farm help
##
gee.ra = geeglm(HelpObserved ~ Distance + RelativeWealthRank + Ego.Size + Alter.Size + Ego.VillageID +
                  HelpReceived,
                id=Ego.HH, family="binomial", data=hh.all, corstr="exchangeable", scale.fix=T)

##
## need
##
gee.need_afb = geeglm(HelpObserved ~ Distance + RelativeWealthRank + Ego.Size + Alter.Size + Ego.VillageID +
                        Ego.Need.AFB + Alter.Need.AFB,
                      id=Ego.HH, family="binomial", data=hh.all, corstr="exchangeable", scale.fix=T)

##
## gifts
##
gee.gift = geeglm(HelpObserved ~ Distance + RelativeWealthRank + Ego.Size + Alter.Size + Ego.VillageID +
                    GiftGiven.ind,
                  id=Ego.HH, family="binomial", data=hh.all, corstr="exchangeable", scale.fix=T)

##
## kin
## (control + partners and kids + r x distance + r x kids)
##
gee.kin = geeglm(HelpObserved ~ Distance + RelativeWealthRank + Ego.Size + Alter.Size + Ego.VillageID +
                   r + r:Distance + AnyWivesInAlter + AnyHusbandsInAlter + ChildrenInAlter + r:ChildrenInAlter,
                 id=Ego.HH, family="binomial", data=hh.all, corstr="exchangeable", scale.fix=T)

##
## kin + need
##
gee.kin_need = geeglm(HelpObserved ~ Distance + RelativeWealthRank + Ego.Size + Alter.Size + Ego.VillageID +
                        r + r:Distance + AnyWivesInAlter + AnyHusbandsInAlter + ChildrenInAlter + r:ChildrenInAlter +
                        Ego.Need.AFB + Alter.Need.AFB,
                      id=Ego.HH, family="binomial", data=hh.all, corstr="exchangeable", scale.fix=T)

##
## kin + reciprocity
##
gee.kin_ra = geeglm(HelpObserved ~ Distance + RelativeWealthRank + Ego.Size + Alter.Size + Ego.VillageID +
                      r + r:Distance + AnyWivesInAlter + AnyHusbandsInAlter + ChildrenInAlter + r:ChildrenInAlter +
                      HelpReceived,
                    id=Ego.HH, family="binomial", data=hh.all, corstr="exchangeable", scale.fix=T)

##
## kin + gifts
##
gee.kin_gift = geeglm(HelpObserved ~ Distance + RelativeWealthRank + Ego.Size + Alter.Size + Ego.VillageID +
                        r + r:Distance + AnyWivesInAlter + AnyHusbandsInAlter + ChildrenInAlter + r:ChildrenInAlter +  # kin
                        GiftGiven.ind,  # gift
                      id=Ego.HH, family="binomial", data=hh.all, corstr="exchangeable", scale.fix=T)

##
## kin + reciprocity + need
##
gee.kin_ra_need = geeglm(HelpObserved ~ Distance + RelativeWealthRank + Ego.Size + Alter.Size + Ego.VillageID +
                           r + r:Distance + AnyWivesInAlter + AnyHusbandsInAlter + ChildrenInAlter + r:ChildrenInAlter +
                           HelpReceived +
                           Ego.Need.AFB + Alter.Need.AFB,
                         id=Ego.HH, family="binomial", data=hh.all, corstr="exchangeable", scale.fix=T)

##
## kin + need + gifts
##
gee.kin_need_gift = geeglm(HelpObserved ~ Distance + RelativeWealthRank + Ego.Size + Alter.Size + Ego.VillageID +
                             r + r:Distance + AnyWivesInAlter + AnyHusbandsInAlter + ChildrenInAlter + r:ChildrenInAlter +  # kin
                             Ego.Need.AFB + Alter.Need.AFB +  # need
                             GiftGiven.ind,  # gift
                           id=Ego.HH, family="binomial", data=hh.all, corstr="exchangeable", scale.fix=T)

##
## kin + need + gifts + recip
##
gee.kin_ra_need_gift = geeglm(HelpObserved ~ Distance + RelativeWealthRank + Ego.Size + Alter.Size + Ego.VillageID +
                                r + r:Distance + AnyWivesInAlter + AnyHusbandsInAlter + ChildrenInAlter + r:ChildrenInAlter +
                                HelpReceived +
                                Ego.Need.AFB + Alter.Need.AFB +
                                GiftGiven.ind,
                              id=Ego.HH, family="binomial", data=hh.all, corstr="exchangeable", scale.fix=T)

##
## kin x need + gifts + recip
##
gee.kin_ra_need_gift_int = geeglm(HelpObserved ~ Distance + RelativeWealthRank + Ego.Size + Alter.Size + Ego.VillageID +
                                    r + r:Distance + AnyWivesInAlter + AnyHusbandsInAlter + ChildrenInAlter + r:ChildrenInAlter +
                                    HelpReceived +
                                    (Ego.Need.AFB + Alter.Need.AFB) * r +
                                    GiftGiven.ind,
                                  id=Ego.HH, family="binomial", data=hh.all, corstr="exchangeable", scale.fix=T)

##
## kin x need + kin x recip + gifts
##
gee.kin_ra_need_gift_int_2 = geeglm(HelpObserved ~ Distance + RelativeWealthRank + Ego.Size + Alter.Size + Ego.VillageID +
                                      r + r:Distance + AnyWivesInAlter + AnyHusbandsInAlter + ChildrenInAlter + r:ChildrenInAlter +
                                      HelpReceived + r:HelpReceived +
                                      (Ego.Need.AFB + Alter.Need.AFB) * r +
                                      GiftGiven.ind,
                                    id=Ego.HH, family="binomial", data=hh.all, corstr="exchangeable", scale.fix=T)

##
## model comparison and averaging
##
best.models = model.sel(gee.null, gee.village, gee.control,
                        gee.kin, gee.ra, gee.need_afb, gee.gift,
                        gee.kin_need, gee.kin_ra, gee.kin_gift,
                        gee.kin_ra_need, gee.kin_need_gift,
                        gee.kin_ra_need_gift, gee.kin_ra_need_gift_int, gee.kin_ra_need_gift_int_2,
                        rank=QIC)
best.models

# average the best (delta QIC < 2) GEE models
# get the top deltaQIC < 2 models and standardize them
top.models   = mget(row.names(subset(best.models, delta<2)))
top.models.z = lapply(top.models, standardize.gee, standardize.y=F, binary.inputs="full")

gee.best   = top.models[[1]]
gee.best.z = top.models.z[[1]]

# average the std and unstd models if there was more than one top models
if (length(top.models) > 1)
{
  gee.best   = model.avg(top.models,   rank=QIC, revised.var=T, fit=T)
  gee.best.z = model.avg(top.models.z, rank=QIC, revised.var=T)
}

save.image(file=file.path(models.dir, "gees - villages.RData"))

max_stars_repo_path: fit gee models.r

max_stars_repo_name: matthewgthomas/mosuo-kinship

max_issues_repo_path: fit gee models.r

max_issues_repo_name: matthewgthomas/mosuo-kinship

max_forks_repo_path: fit gee models.r

max_forks_repo_name: matthewgthomas/mosuo-kinship


text:
context("wald test")

test_that("100 samples, 2 components, same as glm.nb", {
    f1 = runif(100)
    f2 = 1 - f1
    c1 = stats::rnbinom(n=100, size=10, mu=1000)
    c2 = stats::rnbinom(n=100, size=20, mu=10)
    both = round(f1 * c1 + f2 * c2)

    ctl = MASS::glm.nb(both ~ 0 + f1 + f2)
    stat = coef(summary(ctl))[,'z value']

    mf = model.frame(both ~ 0 + f1 + f2)
    terms = attr(mf, "terms")
    x = model.matrix(terms, mf)
    y = model.response(mf)
    res = wald(x, y)

    # expect same coefficients, max 15% diff stat due to diff fitting
    expect_true(all(res$estimate[1:2] - coef(ctl) < 1e-4))
    expect_true(abs(log2(res$statistic[1] / stat[1])) < log2(1.15))
    expect_true(abs(log2(res$statistic[2] / stat[2])) < log2(1.15))
})

max_stars_repo_name: mschubert/boundie

max_issues_repo_name: mschubert/boundie

max_forks_repo_name: mschubert/boundie


text:
# Parallel Growth Power Analysis
# Raymond Viviano
# rayviviano@gmail.com
# February 16th, 2020

suppressPackageStartupMessages(library(lavaan))
suppressPackageStartupMessages(library(mnormt))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(RColorBrewer))

options(warn=-1)

################### CHANGE MODEL AND SCRIPT PARAMETERS HERE ####################
# Measurement model and covariance paths for subjective cognitive ability
scd_intercept_mean_control <- 0
scd_intercept_mean_mci <- 0
scd_slope_mean_control <- 0
scd_slope_mean_mci <- 2 # Vary between 1 and 2
scd_intercept_variance <- 1
scd_slope_variance <- 1
scd_residual_variance <- .5
scd_intercept_slope_cov <- .2
scd_intercept_age_coef_unstd <-  .2
scd_intercept_edu_coef_unstd <- -.2
scd_intercept_sex_coef_unstd <-  .1
scd_intercept_accult_coef_unstd <-  .1
scd_slope_age_coef_unstd <-  .2
scd_slope_edu_coef_unstd <- -.2
scd_slope_sex_coef_unstd <-  .1
scd_slope_accult_coef_unstd <- .1

# Measurement model and covariance paths for objective cognitive performance
obj_intercept_mean_control <- 0
obj_intercept_mean_mci <- 0
obj_slope_mean_control <- 0
# obj_slope_mean_mci varied below in the for loops
obj_intercept_variance <- 1
obj_slope_variance <- 1
obj_residual_variance <- .5
obj_intercept_slope_cov <- .2
obj_intercept_age_coef_unstd <- -.2
obj_intercept_edu_coef_unstd <- .1
obj_intercept_sex_coef_unstd <- .1
obj_intercept_accult_coef_unstd <- .1
obj_slope_age_coef_unstd <- -.2
obj_slope_edu_coef_unstd <- .1
obj_slope_sex_coef_unstd <- .1
obj_slope_accult_coef_unstd <- .1

# Structural model linking growth terms between subjective and objective cog
scd_int_to_obj_int <- -.3
# scd_int_to_obj_slp varied in the for loops
# scd_slp_to_obj_slp varied in the for loops

# Other - General Simulation Parameters
nsimulations <- 5000

#################### END USER DEFINED PARAMETERS SECTION  ######################

# # Variables that the script varies...
# obj_slope_mean_mci; varied between -1 and -2
# scd_int_to_obj_slp; varied between .0 and -1
# scd_slp_to_obj_slp; variedbetween  .0 and -1

sim_df <- data.frame(matrix(ncol = 7, nrow = 0))

colnames(sim_df) <- c("scd_slope_mean_mci", "obj_slope_mean_mci", "scd_int_to_obj_int",
                      "scd_int_to_obj_slp", "scd_slp_to_obj_slp", "rmsea05_si_path05_prop",
                      "rmsea05_ss_path05_prop")

for (obj_slope_mean_mci in c(-1,-2)){
    for (scd_int_to_obj_slp in c(0, -.2, -.4, -.6, -.8, -1)){
        for (scd_slp_to_obj_slp in c(0, -.2, -.4, -.6, -.8, -1)){
            # Measurment model, latent means, vars, and covars for scd latent growth process
            # Reminder, latent slope varies by control/MCI status, so define that later
            growth_model_base <- "i_scd =~ 1*y1 + 1*y2 + 1*y3\n"
            growth_model_base <- paste0(growth_model_base, "s_scd =~ 0*y1 + 1*y2 + 2*y3\n")
            growth_model_base <- paste0(growth_model_base, "i_scd ~~ ", scd_intercept_variance, "*i_scd\n")
            growth_model_base <- paste0(growth_model_base, "s_scd ~~ ", scd_slope_variance, "*s_scd\n")
            growth_model_base <- paste0(growth_model_base, "i_scd ~~ ", scd_intercept_slope_cov, "*s_scd\n")

            # Measurment model, latent means, vars, and covars for objective cognitive
            # performance latent growth process. Reminder, latent slope varies by control
            # or MCI status, so define that later
            growth_model_base <- paste0(growth_model_base, "i_obj =~ 1*x1 + 1*x2 + 1*x3\n")
            growth_model_base <- paste0(growth_model_base, "s_obj =~ 0*x1 + 1*x2 + 2*x3\n")
            growth_model_base <- paste0(growth_model_base, "i_obj ~~ ", obj_intercept_variance, "*i_obj\n")
            growth_model_base <- paste0(growth_model_base, "s_obj ~~ ", obj_slope_variance, "*s_obj\n")
            growth_model_base <- paste0(growth_model_base, "i_obj ~~ ", obj_intercept_slope_cov, "*s_obj\n")

            # Set residual variances for repeated measures of scd
            growth_model_base <- paste0(growth_model_base, "y1 ~~ ", scd_residual_variance, "*y1\n")
            growth_model_base <- paste0(growth_model_base, "y2 ~~ ", scd_residual_variance, "*y2\n")
            growth_model_base <- paste0(growth_model_base, "y3 ~~ ", scd_residual_variance, "*y3\n")

            # Set residual variances for repeated measures of objective performance
            growth_model_base <- paste0(growth_model_base, "x1 ~~ ", obj_residual_variance, "*x1\n")
            growth_model_base <- paste0(growth_model_base, "x2 ~~ ", obj_residual_variance, "*x2\n")
            growth_model_base <- paste0(growth_model_base, "x3 ~~ ", obj_residual_variance, "*x3\n")

            # Age, edu, sex, and acculturation path coefficients to scd latent processes
            growth_model_base <- paste0(growth_model_base, "i_scd ~ ", scd_intercept_age_coef_unstd, "*age\n")
            growth_model_base <- paste0(growth_model_base, "s_scd ~ ", scd_slope_age_coef_unstd, "*age\n")
            growth_model_base <- paste0(growth_model_base, "i_scd ~ ", scd_intercept_edu_coef_unstd, "*edu\n")
            growth_model_base <- paste0(growth_model_base, "s_scd ~ ", scd_slope_edu_coef_unstd, "*edu\n")
            growth_model_base <- paste0(growth_model_base, "i_scd ~ ", scd_intercept_sex_coef_unstd, "*sex\n")
            growth_model_base <- paste0(growth_model_base, "s_scd ~ ", scd_slope_sex_coef_unstd, "*sex\n")
            growth_model_base <- paste0(growth_model_base, "i_scd ~ ", scd_intercept_accult_coef_unstd, "*acculturation\n")
            growth_model_base <- paste0(growth_model_base, "s_scd ~ ", scd_slope_accult_coef_unstd, "*acculturation\n")

            # Age, edu, sex, and acculturation path coefficients to objective
            # cognitive performance latent processes
            growth_model_base <- paste0(growth_model_base, "i_obj ~ ", obj_intercept_age_coef_unstd, "*age\n")
            growth_model_base <- paste0(growth_model_base, "s_obj ~ ", obj_slope_age_coef_unstd, "*age\n")
            growth_model_base <- paste0(growth_model_base, "i_obj ~ ", obj_intercept_edu_coef_unstd, "*edu\n")
            growth_model_base <- paste0(growth_model_base, "s_obj ~ ", obj_slope_edu_coef_unstd, "*edu\n")
            growth_model_base <- paste0(growth_model_base, "i_obj ~ ", obj_intercept_sex_coef_unstd, "*sex\n")
            growth_model_base <- paste0(growth_model_base, "s_obj ~ ", obj_slope_sex_coef_unstd, "*sex\n")
            growth_model_base <- paste0(growth_model_base, "i_obj ~ ", obj_intercept_accult_coef_unstd, "*acculturation\n")
            growth_model_base <- paste0(growth_model_base, "s_obj ~ ", obj_slope_accult_coef_unstd, "*acculturation\n")

            # Create group specific models from the base model
            growth_model_con <- growth_model_base
            growth_model_mci <- growth_model_base

            # Set group specific paramters for different group models
            growth_model_con <- paste0(growth_model_con, "i_scd ~ ", scd_intercept_mean_control, "*1\n")
            growth_model_mci <- paste0(growth_model_mci, "i_scd ~ ", scd_intercept_mean_mci, "*1\n")
            growth_model_con <- paste0(growth_model_con, "s_scd ~ ", scd_slope_mean_control, "*1\n")
            growth_model_mci <- paste0(growth_model_mci, "s_scd ~ ", scd_slope_mean_mci, "*1\n")
            growth_model_con <- paste0(growth_model_con, "i_obj ~ ", obj_intercept_mean_control, "*1\n")
            growth_model_mci <- paste0(growth_model_mci, "i_obj ~ ", obj_intercept_mean_mci, "*1\n")
            growth_model_con <- paste0(growth_model_con, "s_obj ~ ", obj_slope_mean_control, "*1\n")
            growth_model_mci <- paste0(growth_model_mci, "s_obj ~ ", obj_slope_mean_mci, "*1\n")

            # Define structural model. Add paths between growth terms for con and mci models
            growth_model_con <- paste0(growth_model_con, "i_obj ~ ", scd_int_to_obj_int, "*i_scd\n")
            growth_model_mci <- paste0(growth_model_mci, "i_obj ~ ", scd_int_to_obj_int, "*i_scd\n")
            growth_model_con <- paste0(growth_model_con, "s_obj ~ ", scd_int_to_obj_slp, "*i_scd\n")
            growth_model_mci <- paste0(growth_model_mci, "s_obj ~ ", scd_int_to_obj_slp, "*i_scd\n")
            growth_model_con <- paste0(growth_model_con, "s_obj ~ ", scd_slp_to_obj_slp, "*s_scd\n")
            growth_model_mci <- paste0(growth_model_mci, "s_obj ~ ", scd_slp_to_obj_slp, "*s_scd\n")

            # Test Model - Correctly Specified
            m1 <- "
            i_scd =~ 1*y1 + 1*y2 + 1*y3
            s_scd =~ 0*y1 + 1*y2 + 2*y3
            i_scd ~~ i_scd + s_scd
            s_scd ~~ s_scd

            i_obj =~ 1*x1 + 1*x2 + 1*x3
            s_obj =~ 0*x1 + 1*x2 + 2*x3
            i_obj ~~ i_obj + s_obj
            s_obj ~~ s_obj

            i_scd ~ mci + age + edu + sex + acculturation
            s_scd ~ mci + age + edu + sex + acculturation
            i_obj ~ mci + age + edu + sex + acculturation + i_scd
            s_obj ~ mci + age + edu + sex + acculturation + i_scd + s_scd
            "

            # Run the simulations, SI <- Slope ON Intercept, SS <- Slope ON Slope
            rmsea05_si_path05 <- 0
            rmsea05_ss_path05 <- 0

            # Keep track of successful simulations without convergence issues
            ngoodsims <- 0

            for(k in 1:nsimulations){
                # Simulate 100 MCI and 100 Control Observations
                df_con <- lavaan::simulateData(model=growth_model_con, model.type="growth",
                                                sample.nobs=100)

                df_mci <- lavaan::simulateData(model=growth_model_mci, model.type="growth",
                                                sample.nobs=100)

                # Add MCI columns to both data frames
                df_con[,"mci"] <- 0
                df_mci[,"mci"] <- 1

                # Combine Control and MCI datasets
                df <- rbind(df_con, df_mci)

                # Test the correctly specified model
                e <- try(fit <- lavaan::growth(m1, data=df, missing='fiml'))

                # Skip this simulation if model did not converge
                if(inherits(e, "try-error")) {
                    fit <- NULL
                }

                # Check fit errors/model convergence errors
                e2 <- try(fitStatsDummy <- fitmeasures(fit))
                if(inherits(e2, "try-error")) {
                    fit <- NULL
                }

                # If the model ran successfully, extract fit stats and parameter estimates
                if (!is.null(fit)) {
                    # Increment number of successful simulations
                    ngoodsims <- ngoodsims + 1

                    # Get RMSEA, s_obj~i_scd and s_obj~s_scd pvalues
                    rmsea <- as.double(fitMeasures(fit, c('rmsea'))[1])
                    si_pval <- as.double(parameterEstimates(fit, standardized = FALSE)[40, 'pvalue'])
                    ss_pval <- as.double(parameterEstimates(fit, standardized = FALSE)[41, 'pvalue'])

                    if(rmsea <= .05){
                        if(si_pval <= .05){
                            rmsea05_si_path05 <- rmsea05_si_path05 + 1
                        } else if(ss_pval <= .05){
                        rmsea05_ss_path05 <- rmsea05_ss_path05 + 1
                        }
                    }
                }
            }

            rmsea05_si_path05_prop <- rmsea05_si_path05/ngoodsims
            rmsea05_ss_path05_prop <- rmsea05_ss_path05/ngoodsims

            # Create dataframe to store simulation results
            temp <- data.frame(scd_slope_mean_mci, obj_slope_mean_mci, scd_int_to_obj_int,
                                    scd_int_to_obj_slp, scd_slp_to_obj_slp, rmsea05_si_path05_prop,
                                    rmsea05_ss_path05_prop)

            colnames(temp) <- c("scd_slope_mean_mci", "obj_slope_mean_mci", "scd_int_to_obj_int",
                                "scd_int_to_obj_slp", "scd_slp_to_obj_slp", "rmsea05_si_path05_prop",
                                "rmsea05_ss_path05_prop")

            sim_df <- rbind(sim_df, temp)

        }
    }
}

foo <- as.character(Sys.time())

print(foo)
print(paste0("./sim_parallel_growth-", as.character(Sys.time()),".Rda"))

save(sim_df, file=paste0("./sim_parallel_growth-", as.character(Sys.time()),".Rda"))
write.csv(sim_df, paste0("./sim_parallel_growth-", as.character(Sys.time()),".csv"))

max_stars_repo_name: rviviano/data-tools

max_issues_repo_name: rviviano/data-tools

max_forks_repo_name: rviviano/data-tools