simNormMat <- function(n1, n2, t, alpha = 0.05, nsims = 1000)
{
nobs <- (n1+n2)*nsims
simMat <- matrix(NA,nrow=nobs,ncol=nsims)
for (i in 1:nsims)
{
samp1 <- rnorm(n1)
samp2 <- rnorm(n2)
m <- min(n1,n2)
idx1 <- sample(1:m,t)
idx2 <- sample(1:m,t)
samp1[idx1] <- samp2[idx2]
simMat[,i] <- c(samp1,samp2)
}
return(simMat)
}
simNormMat(4,4,1,nsims=10)
simNormMat <- function(n1, n2, t, alpha = 0.05, nsims = 1000)
{
totalN <- n1 + n2
nobs <- totalN*nsims
simMat <- matrix(NA,nrow=totalN,ncol=nsims)
for (i in 1:nsims)
{
samp1 <- rnorm(n1)
samp2 <- rnorm(n2)
m <- min(n1,n2)
idx1 <- sample(1:m,t)
idx2 <- sample(1:m,t)
samp1[idx1] <- samp2[idx2]
simMat[,i] <- c(samp1,samp2)
}
return(simMat)
}
simNormMat(4,4,1,nsims=10)
funToRepeat <- function(dims) {
err <- runif(prod(dims))
matrix(err,nrow=dims[1],ncol=dims[2])
}
library(dplyr)
ll <- alply(cbind(rep(10,10),rep(20,10)),1,funToRepeat)
install.packages("plyr")
library(plyr)
ll <- alply(cbind(rep(10,10),rep(20,10)),1,funToRepeat)
ll
funToRepeat(2,2)
funToRepeat(2)
funToRepeat(c(2,2)
)
funToRepeat
help(alply)
wmwSimNormalties <- function(n1, n2, t, ties.method, alpha = 0.05, nsims = 1000)
{	totalN <- n1+n2
nobs <- totalN*nsims
simMat <- simNormMat(n1,n2,t,nsims)
dat <- data.frame(Group=gl(2,n1,length=totalN),simMat)
WMWpvals <- sapply(dat[,2:ncol(dat)],function(xx) pvalue(wilcox_test(xx~dat$Group, ties.method=ties.method)))
Tpvals <- sapply(dat[,2:ncol(dat)],function(xx) t.test(xx~dat$Group)$p.value)
WMWerror <- length(WMWpvals[WMWpvals < alpha])/nsims
Terror <- length(Tpvals[Tpvals < alpha])/nsims
return(data.frame(error.WMW=WMWerror, error.t = Terror))
}
wmwSimNormalties(10,10,2,"mid-ranks",nsims=10)
simNormMat <- function(n1, n2, t, alpha = 0.05, nsims = 1000)
{
totalN <- n1 + n2
simMat <- matrix(NA,nrow=totalN,ncol=nsims)
for (i in 1:nsims)
{
samp1 <- rnorm(n1)
samp2 <- rnorm(n2)
m <- min(n1,n2)
idx1 <- sample(1:m,t)
idx2 <- sample(1:m,t)
samp1[idx1] <- samp2[idx2]
simMat[,i] <- c(samp1,samp2)
}
return(simMat)
}
wmwSimNormalties(10,10,2,"mid-ranks",nsims=10)
wmwSimNormalties <- function(n1, n2, t, ties.method, alpha = 0.05, nsims = 1000)
{	totalN <- n1+n2
nobs <- totalN*nsims
simMat <- simNormMat(n1,n2,t,nsims)
dat <- data.frame(Group=gl(2,n1,length=totalN),simMat)
WMWpvals <- sapply(dat[,2:ncol(dat)],function(xx) pvalue(wilcox_test(xx~dat$Group,ties.method=ties.method)))
Tpvals <- sapply(dat[,2:ncol(dat)],function(xx) t.test(xx~dat$Group)$p.value)
WMWerror <- length(WMWpvals[WMWpvals < alpha])/nsims
Terror <- length(Tpvals[Tpvals < alpha])/nsims
return(data.frame(error.WMW=WMWerror, error.t = Terror))
}
library(coin)
library(VGAM)
wmwSimNormalties(10,10,2,"mid-ranks",nsims=10)
help("data.frame")
wmwSimNormalties <- function(n1, n2, t, ties.method, alpha = 0.05, nsims = 1000)
{	totalN <- n1+n2
nobs <- totalN*nsims
simMat <- simNormMat(n1,n2,t,nsims)
dat <- data.frame(Group=gl(2,n1,length=totalN),simMat)
WMWpvals <- sapply(dat[,2:ncol(dat)],function(xx) pvalue(wilcox_test(xx~dat$Group,ties.method=ties.method)))
Tpvals <- sapply(dat[,2:ncol(dat)],function(xx) t.test(xx~dat$Group)$p.value)
WMWerror <- length(WMWpvals[WMWpvals < alpha])/nsims
Terror <- length(Tpvals[Tpvals < alpha])/nsims
return(data.frame(test=c("WMW","Ttest"),error=c(WMWerror, Terror))
}
wmwSimNormalties <- function(n1, n2, t, ties.method, alpha = 0.05, nsims = 1000)
{	totalN <- n1+n2
nobs <- totalN*nsims
simMat <- simNormMat(n1,n2,t,nsims)
dat <- data.frame(Group=gl(2,n1,length=totalN),simMat)
WMWpvals <- sapply(dat[,2:ncol(dat)],function(xx) pvalue(wilcox_test(xx~dat$Group,ties.method=ties.method)))
Tpvals <- sapply(dat[,2:ncol(dat)],function(xx) t.test(xx~dat$Group)$p.value)
WMWerror <- length(WMWpvals[WMWpvals < alpha])/nsims
Terror <- length(Tpvals[Tpvals < alpha])/nsims
return(data.frame(test=c("WMW","Ttest"),error=c(WMWerror, Terror)))
}
wmwSimNormalties(10,10,2,"mid-ranks",nsims=10)
help("seq_len")
lapply(seq_len(2), wmwSimNormalties(10,10,2,"mid-ranks",nsims=10))
gen_mat <- function(x) matrix(c(1, 1, 1, x + rnorm(1)), nrow = 2)
gen_mat(1)
gen_mat(3)
n <- 10
lapply(seq_len(n), gen_mat(1))
help(lapply(list, function))
help(lapply
)
sample(letters(1:15))
sample(letters[1:4],15)
sample(letters[1:4],15,replace=TRUE)
sample(letters[1:4],15,replace=TRUE)
10.9 +0.23*73
29 - (10.9 +0.23*73)
27.69 - 29
10.9 +0.23*70
library(MASS)
data(birthwt)
str(birthwt)
library(fpp)
install.packages("fpp2")
install.packages("fpp2")
library(fpp2)
sessionInfo()
library(fpp2)
data(ausair)
air <- window(ausair,start=1990, end=2004)
fit1 <- holt(air, alpha=0.8, beta=0.2, initial="simple", h=5)
accuracy(fit1,ausair)
fit2 <- holt(air, alpha=0.8, beta=0.2, initial="simple", exponential=TRUE, h=5)
accuracy(fit2,ausair)
fit3 <- holt(air, alpha=0.8, beta=0.2, damped=TRUE, initial="simple", h=5)
fit3 <- holt(air, alpha=0.8, beta=0.2, damped=TRUE, h=5)
accuracy(fit3,ausair)
help("accuracy")
library(fpp2)
data(usdeaths)
data("sunspotarea")
autoplot(usdeaths)
ggseasonplot(usdeaths)
ggsubseriesplot(usdeaths)
autoplot(sunspotarea)
ggseasonplot(sunspotarea)
ggsubseriesplot(sunspotarea)
sample(1:19,19)
setwd()
getwd()
setwd("./Dropbox/2019 Spring/")
list.files()
load("Combined_Jan14_SimResults.Rdata")
setting <- expand.grid(list("p" = c(10 , 50 , 100 , 200) ,
"amt_ME" = c(.1 , .25 , .5) ,
"nMCMC" = c(20 , 50 , 100) ,
"graphType" = c("hub" , "random" , "band" , "cluster")))
setting[80,]
setting
simulationOutput[[55]]
simulationOutput[[56]]
simulationOutput[[127]]
simulationOutput[[131]]
setwd("~/Dropbox/2019 Spring/FlowThroughCentrality/Simulations")
knitr::opts_chunk$set(echo = TRUE)
setwd("Nodes50Deg6Pert6")
setwd("~/Dropbox/2019 Spring/FlowThroughCentrality/Simulations")
getwd()
setwd("Nodes50Deg6Pert3")
setwd("Nodes50Deg6Pert3")
getwd()
setwd("/Nodes50Deg6Pert3")
setwd("~/Nodes50Deg6Pert3")
setwd("./Nodes50Deg6Pert3")
varNames <- c("Network",paste0(c("Cent","Tier"), rep(1:50,each=2)))
dataBW <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert03_Betweenness.csv",skip=5,nrows=100,header=FALSE)
colnames(dataBW) <- varNames
dataBW$Measure <- rep("BW",100)
dataCL <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert03_Closeness.csv",skip=5,nrows=100,header=FALSE)
colnames(dataCL) <- varNames
dataCL$Measure <- rep("CL",100)
dataFB <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert03_FlowBetween.csv",skip=5,nrows=100,header=FALSE)
colnames(dataFB) <- varNames
dataFB$Measure <- rep("FB",100)
dataFT <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert03_Flowthrough.csv",skip=5,nrows=100,header=FALSE)
colnames(dataFT) <- varNames
dataFT$Measure <- rep("FT",100)
dataSB <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert03_StableBetween.csv",skip=5,nrows=100,header=FALSE)
colnames(dataSB) <- varNames
dataSB$Measure <- rep("SB",100)
## Create tibbles with first six columns of each file
dataBWfirst6 <- dplyr::select(dataBW,c(Network:Tier6,Measure))
dataCLfirst6 <- dplyr::select(dataCL,c(Network:Tier6,Measure))
dataFBfirst6 <- dplyr::select(dataFB,c(Network:Tier6,Measure))
dataFTfirst6 <- dplyr::select(dataFT,c(Network:Tier6,Measure))
dataSBfirst6 <- dplyr::select(dataSB,c(Network:Tier6,Measure))
dataAllFirst6 <- bind_rows(dataBWfirst6,dataCLfirst6,dataFBfirst6,dataFTfirst6, dataSBfirst6)
library(coin) # load the package in R
suppressPackageStartupMessages(library(tidyverse)) # for tidying the data
library(xtable) # for obtaining tables in LaTex format
varNames <- c("Network",paste0(c("Cent","Tier"), rep(1:50,each=2)))
dataBW <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert03_Betweenness.csv",skip=5,nrows=100,header=FALSE)
colnames(dataBW) <- varNames
dataBW$Measure <- rep("BW",100)
dataCL <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert03_Closeness.csv",skip=5,nrows=100,header=FALSE)
colnames(dataCL) <- varNames
dataCL$Measure <- rep("CL",100)
dataFB <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert03_FlowBetween.csv",skip=5,nrows=100,header=FALSE)
colnames(dataFB) <- varNames
dataFB$Measure <- rep("FB",100)
dataFT <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert03_Flowthrough.csv",skip=5,nrows=100,header=FALSE)
colnames(dataFT) <- varNames
dataFT$Measure <- rep("FT",100)
dataSB <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert03_StableBetween.csv",skip=5,nrows=100,header=FALSE)
colnames(dataSB) <- varNames
dataSB$Measure <- rep("SB",100)
## Create tibbles with first six columns of each file
dataBWfirst6 <- dplyr::select(dataBW,c(Network:Tier6,Measure))
dataCLfirst6 <- dplyr::select(dataCL,c(Network:Tier6,Measure))
dataFBfirst6 <- dplyr::select(dataFB,c(Network:Tier6,Measure))
dataFTfirst6 <- dplyr::select(dataFT,c(Network:Tier6,Measure))
dataSBfirst6 <- dplyr::select(dataSB,c(Network:Tier6,Measure))
dataAllFirst6 <- bind_rows(dataBWfirst6,dataCLfirst6,dataFBfirst6,dataFTfirst6, dataSBfirst6)
dataAllFirst6$Measure <- as.factor(dataAllFirst6$Measure)
tier6n50deg6pert3a <- gather(dataAllFirst6,key=CentLevel,value=Centrality,Cent1,Cent2,Cent3,Cent4,Cent5,Cent6,factor_key=TRUE)
tier6n50deg6pert3b <- gather(dataAllFirst6,key=Level,value=Tier, Tier1,Tier2,Tier3,Tier4,Tier5,Tier6, factor_key=TRUE)
tier6n50deg6pert3 <-merge(tier6n50deg6pert3a,tier6n50deg6pert3b,by=c("Network","Measure"))
tier6n50deg6pert3 <- dplyr::select(tier6n50deg6pert3,-c(Tier1:Tier3,Cent1:Cent3))
tier6n50deg6pert3 <- separate(tier6n50deg6pert3,Network,into=c("Trash","Sample"),sep=" ")
tier6n50deg6pert3$Sample <- as.integer(tier6n50deg6pert3$Sample)
tier6n50deg6pert3 <- tier6n50deg6pert3[,-1]
setwd("./Nodes50Deg6Pert6")
setwd("../Nodes50Deg6Pert6")
getwd()
setwd("../Nodes50Deg6Pert6")
varNames <- c("Network",paste0(c("Cent","Tier"), rep(1:50,each=2)))
dataBW <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert06_Betweenness.csv",skip=5,nrows=100,header=FALSE)
getwd()
dataBW <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert06_Betweenness.csv",skip=5,nrows=100,header=FALSE)
colnames(dataBW) <- varNames
dataBW$Measure <- rep("BW",100)
dataCL <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert06_Closeness.csv",skip=5,nrows=100,header=FALSE)
dataCL <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert06_Closeness.csv",skip=5,nrows=100,header=FALSE)
colnames(dataCL) <- varNames
dataCL$Measure <- rep("CL",100)
dataFB <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert06_FlowBetween.csv",skip=5,nrows=100,header=FALSE)
dataFB <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert06_FlowBetween.csv",skip=5,nrows=100,header=FALSE)
colnames(dataFB) <- varNames
dataFB$Measure <- rep("FB",100)
dataFT <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert06_Flowthrough.csv",skip=5,nrows=100,header=FALSE)
dataFT <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert06_Flowthrough.csv",skip=5,nrows=100,header=FALSE)
colnames(dataFT) <- varNames
dataFT$Measure <- rep("FT",100)
dataSB <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert06_StableBetween.csv",skip=5,nrows=100,header=FALSE)
dataSB <- read.csv(file="Nodes50_AvgDeg06_UW_01_Pert06_StableBetween.csv",skip=5,nrows=100,header=FALSE)
colnames(dataSB) <- varNames
dataSB$Measure <- rep("SB",100)
dataBWfirst6 <- dplyr::select(dataBW,c(Network:Tier6,Measure))
dataCLfirst6 <- dplyr::select(dataCL,c(Network:Tier6,Measure))
dataFBfirst6 <- dplyr::select(dataFB,c(Network:Tier6,Measure))
dataFTfirst6 <- dplyr::select(dataFT,c(Network:Tier6,Measure))
dataSBfirst6 <- dplyr::select(dataSB,c(Network:Tier6,Measure))
dataAllFirst6 <- bind_rows(dataBWfirst6,dataCLfirst6,dataFBfirst6,dataFTfirst6, dataSBfirst6)
dataAllFirst6$Measure <- as.factor(dataAllFirst6$Measure)
tier6n50deg6pert6a <- gather(dataAllFirst6,key=CentLevel,value=Centrality,Cent1,Cent2,Cent3,Cent4,Cent5,Cent6,factor_key=TRUE)
tier6n50deg6pert6b <- gather(dataAllFirst6,key=Level,value=Tier, Tier1,Tier2,Tier3,Tier4,Tier5,Tier6, factor_key=TRUE)
tier6n50deg6pert6 <-merge(tier6n50deg6pert6a,tier6n50deg6pert6b,by=c("Network","Measure"))
tier6n50deg6pert6 <- dplyr::select(tier6n50deg6pert6,-c(Tier1:Tier3,Cent1:Cent3))
tier6n50deg6pert6 <- separate(tier6n50deg6pert6,Network,into=c("Trash","Sample"),sep=" ")
tier6n50deg6pert6$Sample <- as.integer(tier6n50deg6pert6$Sample)
tier6n50deg6pert6 <- tier6n50deg6pert6[,-1]
tier6n50deg6pert3$Perturbation <- "3%"
tier6n50deg6pert6$Perturbation <- "6%"
n50deg6all <- cbind(tier6n50deg6pert3,tier6n50deg6pert6)
n50deg6all$Perturbation <- as.factor(n50deg6all$Perturbation)
names(n50deg6all)
n50deg6all <- bind_rows(tier6n50deg6pert3,tier6n50deg6pert6)
names(n50deg6all)
names(dataBW)
names(dataBWfirst6)
names(dataAllFirst6)
names(tier6n50deg6pert3a)
names(tier6n50deg6pert3b)
names(tier6n50deg6pert3)
tier6n50deg6pert3 <-merge(tier6n50deg6pert3a,tier6n50deg6pert3b,by=c("Network","Measure"))
names(tier6n50deg6pert3)
summary(tier6n50deg6pert3)
tier6n50deg6pert3$Level
tier6n50deg6pert3$CentLevel
n50deg6all %>% group_by(Perturbation,CentLevel,Measure) %>% summarise(mean=mean(Centrality),sd=sd(Centrality),cv=100*(sd(Centrality)/mean(Centrality)))
n50deg6all %>% group_by(Perturbation,CentLevel,Measure) %>% summarise(mean=round(mean(Centrality),3),sd=round(sd(Centrality),3),cv=round(100*(sd(Centrality)/mean(Centrality)),3))
help(spread)
spread(n50deg6all,key=Perturbation,value)
spread(n50deg6all,key=Perturbation,value=mean)
n50deg6all %>% group_by(Perturbation,CentLevel,Measure) %>% summarise(Mean=round(mean(Centrality),3),SD=round(sd(Centrality),3),CV=round(100*(sd(Centrality)/mean(Centrality)),3))
spread(n50deg6all,key=Perturbation,value=Mean)
n50deg6all %>% group_by(Perturbation,CentLevel,Measure) %>% summarise(Mean=round(mean(Centrality),3),SD=round(sd(Centrality),3),CV=round(100*(sd(Centrality)/mean(Centrality)),3))
spread(n50deg6all,key=Perturbation,value=Mean)
names(n50deg6all)
centSummary <- n50deg6all %>% group_by(Perturbation,CentLevel,Measure) %>% summarise(Mean=round(mean(Centrality),3),SD=round(sd(Centrality),3),CV=round(100*(sd(Centrality)/mean(Centrality)),3))
spread(centSummary,key=Perturbation,value=Mean)
help("unite")
centGather <- centSummary %>% gather(Mean,SD,CV)
names(centGaether)
names(centGather)
centGather <- centSummary %>% gather(Mean,SD)
centGather <- centSummary %>% gather(Mean,SD,CV)
spread(centGather,Perturbation,Mean)
head(centGather)
centGather <- centSummary %>% gather(Mean,CV)
centGather <- centSummary %>% gather(Perturbation,Mean,SD,CV)
head(centGather)
