########################################
#
#   MoNAn workshop Oct 2023, Zurich
#   Script 3
#
########################################

# data preparation from script 1

N_ind <- nrow(mobilityEdgelist)
N_org <- length(unique(as.numeric(mobilityEdgelist)))

transfers <- createEdgelist(mobilityEdgelist,
                            nodeSet = c("organisations", "organisations", "people")
)
people <- createNodeSet(N_ind)
organisations <- createNodeSet(N_org)
sameRegion <- outer(orgRegion, orgRegion, "==") * 1
sameRegion <- createNetwork(sameRegion,
                            nodeSet = c("organisations", "organisations")
)
region <- createNodeVariable(orgRegion, nodeSet = "organisations")
size <- createNodeVariable(orgSize, nodeSet = "organisations", addSim = TRUE)
sex <- createNodeVariable(indSex, nodeSet = "people")

myState <- createProcessState(
  list(
    transfers = transfers,
    people = people,
    organisations = organisations,
    sameRegion = sameRegion,
    region = region,
    size = size,
    sex = sex), 
  dependentVariable = "transfers")

myCache <- createWeightedCache(myState, resourceCovariates = c("sex"))


##### score-type tests #####
# full model from script 1
#
# myResDN
# Results
# Effects   Estimates StandardErrors  Convergence
# 1                                   loops  2.59297741     0.17458651  0.071528769
# 2                         reciprocity_min  0.81086410     0.17908430 -0.044339091
# 3             dyadic_covariate sameRegion  1.68707009     0.11428882 -0.090555885
# 4                    alter_covariate size  0.03749522     0.02341118  0.030002067
# 5 resource_covar_to_node_covar region sex -0.63677043     0.17178533 -0.008664495
# 6                loops_resource_covar sex -0.36464266     0.21616269 -0.004067316
# 
# 
# compared to script 1, we exclude "loops_resource_covar sex" (n.s.) and
# "dyadic_covariate sameRegion" (p-val < 0.001)

myEffects_red <- createEffectsObject(
  list(
    list("loops"),
    list("reciprocity_min"),
    list("alter_covariate", attribute.index = "size"),
    list("resource_covar_to_node_covar",
         attribute.index = "region",
         resource.attribute.index = "sex"
    )
  ), 
  checkProcessState = myState  
)

myAlg_red <- createAlgorithm(myState, myEffects_red, multinomialProposal = FALSE,
                         initGain = 0.8)

myResDN_red <- estimateMobilityNetwork(
  myState, myCache, myEffects_red, myAlg_red,
  parallel = TRUE, cpus = 4,
  verbose = TRUE,
  # The following is very important!!!
  returnDeps = TRUE,
  fish = FALSE
)

myResDN_red

# compare to the full model from script 1
#
# myResDN
# Results
# Effects   Estimates StandardErrors  Convergence
# 1                                   loops  2.59297741     0.17458651  0.071528769
# 2                         reciprocity_min  0.81086410     0.17908430 -0.044339091
# 3             dyadic_covariate sameRegion  1.68707009     0.11428882 -0.090555885
# 4                    alter_covariate size  0.03749522     0.02341118  0.030002067
# 5 resource_covar_to_node_covar region sex -0.63677043     0.17178533 -0.008664495
# 6                loops_resource_covar sex -0.36464266     0.21616269 -0.004067316


### illustration of score tests ###
# extract first simulation from phase 3
myResDN_red$deps[[1]]

# where is the network representation stored?
sim_net <- myResDN_red$deps[[1]]$cache$transfers$valuedNetwork

# calculate same region statistic for this by hand
sum(sameRegion$data * sim_net)

# compare to observation
sum(sameRegion$data * myCache$transfers$valuedNetwork)

# there seems to be quite a big difference...

### now automate this over all simulations
candidate_effects <- createEffectsObject(
  list(
    list("dyadic_covariate", attribute.index = "sameRegion"),
    list("loops_resource_covar", resource.attribute.index = "sex")
  )
)

myTest <- scoreTest(myResDN_red, candidate_effects)

myTest


### other candidate effects

candidate_effects_2 <- createEffectsObject(
  list(
    list("concentration_basic"),
    list("concentration_GW"),
    list("transitivity_min")
  )
)

myTest_2 <- scoreTest(myResDN_red, candidate_effects_2)

myTest_2

# note that this does not warn against concentration_basic as degenerate!


##### degree controls #####
### option 1: demand as a node covariate

orgIndegree <- log(table(mobilityEdgelist[,2])) - mean(log(table(mobilityEdgelist[,2])))
orgIndegree <- table(mobilityEdgelist[,2]) - mean(table(mobilityEdgelist[,2]))
orgDemand <- createNodeVariable(orgIndegree, nodeSet = "organisations")

myState_degree_1 <- createProcessState(
  list(
    transfers = transfers,
    people = people,
    organisations = organisations,
    sameRegion = sameRegion,
    # here
    orgDemand = orgDemand), 
  dependentVariable = "transfers")

myCache_degree_1 <- createWeightedCache(myState_degree_1)

myEffects_degree_1 <- createEffectsObject(
  list(
    list("alter_covariate", attribute.index = "orgDemand")
  ), 
  checkProcessState = myState  
)

myAlg_degree_1 <- createAlgorithm(myState, myEffects_degree_1, multinomialProposal = FALSE,
                                  initGain = 0.8)

myResDN_degree_1 <- estimateMobilityNetwork(
  myState_degree_1, myCache_degree_1, myEffects_degree_1, myAlg_degree_1,
  parallel = TRUE, cpus = 4,
  verbose = TRUE,
  returnDeps = TRUE,
  fish = FALSE
)

myResDN_degree_1

# note that the parameter is very close to 1 (because we used the log)

myEffects_degree_2 <- createEffectsObject(
  list(
    list("loops"),
    list("reciprocity_min"),
    list("dyadic_covariate", attribute.index = "sameRegion"),
    list("alter_covariate", attribute.index = "orgDemand")
  ), 
  checkProcessState = myState  
)


myAlg_degree_2 <- createAlgorithm(myState_degree_1, myEffects_degree_2, multinomialProposal = FALSE,
                                  initGain = 0.8)

myResDN_degree_2 <- estimateMobilityNetwork(
  myState_degree_1, myCache_degree_1, myEffects_degree_2, myAlg_degree_2,
  parallel = TRUE, cpus = 4,
  verbose = TRUE,
  returnDeps = TRUE,
  fish = FALSE
)

myResDN_degree_2

# note that given the other parameters, the estimate for the degree control is not 1 anymore!


### option 2: fixed effects

# create a dummy variable for org 2 to 17
createDummy <- function(length, position){
  vec <- rep(0, length)
  vec[position] <- 1
  vec
}

dummies <- lapply(2:17, function(i) createDummy(17, i))


# create a state including all dummies
# note that the node covariate are created inside the createProcessState function!
myState_fe <- createProcessState(
  list(
    transfers = transfers,
    people = people,
    organisations = organisations,
    sameRegion = sameRegion,
    # here
    dummy2 = createNodeVariable(as.numeric(dummies[[1]]), nodeSet = 'organisations'), 
    dummy3 = createNodeVariable(as.numeric(dummies[[2]]), nodeSet = 'organisations'), 
    dummy4 = createNodeVariable(as.numeric(dummies[[3]]), nodeSet = 'organisations'), 
    dummy5 = createNodeVariable(as.numeric(dummies[[4]]), nodeSet = 'organisations'), 
    dummy6 = createNodeVariable(as.numeric(dummies[[5]]), nodeSet = 'organisations'), 
    dummy7 = createNodeVariable(as.numeric(dummies[[6]]), nodeSet = 'organisations'), 
    dummy8 = createNodeVariable(as.numeric(dummies[[7]]), nodeSet = 'organisations'), 
    dummy9 = createNodeVariable(as.numeric(dummies[[8]]), nodeSet = 'organisations'), 
    dummy10 = createNodeVariable(as.numeric(dummies[[9]]), nodeSet = 'organisations'), 
    dummy11 = createNodeVariable(as.numeric(dummies[[10]]), nodeSet = 'organisations'), 
    dummy12 = createNodeVariable(as.numeric(dummies[[11]]), nodeSet = 'organisations'), 
    dummy13 = createNodeVariable(as.numeric(dummies[[12]]), nodeSet = 'organisations'), 
    dummy14 = createNodeVariable(as.numeric(dummies[[13]]), nodeSet = 'organisations'), 
    dummy15 = createNodeVariable(as.numeric(dummies[[14]]), nodeSet = 'organisations'), 
    dummy16 = createNodeVariable(as.numeric(dummies[[15]]), nodeSet = 'organisations'), 
    dummy17 = createNodeVariable(as.numeric(dummies[[16]]), nodeSet = 'organisations') 
  ), 
  dependentVariable = "transfers")

myCache_fe <- createWeightedCache(myState_fe)


# include all dummies in the effects object
myEffects_fe <- createEffectsObject(
  list(
    list("loops"),
    list("reciprocity_min"),
    list("dyadic_covariate", attribute.index = "sameRegion"),
    
    list('alter_covariate', attribute.index = 'dummy2'),
    list('alter_covariate', attribute.index = 'dummy3'),
    list('alter_covariate', attribute.index = 'dummy4'),
    list('alter_covariate', attribute.index = 'dummy5'),
    list('alter_covariate', attribute.index = 'dummy6'),
    list('alter_covariate', attribute.index = 'dummy7'),
    list('alter_covariate', attribute.index = 'dummy8'),
    list('alter_covariate', attribute.index = 'dummy9'),
    list('alter_covariate', attribute.index = 'dummy10'),
    list('alter_covariate', attribute.index = 'dummy11'),
    list('alter_covariate', attribute.index = 'dummy12'),
    list('alter_covariate', attribute.index = 'dummy13'),
    list('alter_covariate', attribute.index = 'dummy14'),
    list('alter_covariate', attribute.index = 'dummy15'),
    list('alter_covariate', attribute.index = 'dummy16'),
    list('alter_covariate', attribute.index = 'dummy17')
  ), 
  checkProcessState = myState_fe 
)

myAlg_fe <- createAlgorithm(myState_fe, myEffects_fe, multinomialProposal = FALSE,
                                  initGain = 0.8)

# careful, this now takes a lot longer!
myResDN_fe <- estimateMobilityNetwork(
  myState_fe, myCache_fe, myEffects_fe, myAlg_fe,
  parallel = TRUE, cpus = 8,
  verbose = TRUE,
  returnDeps = FALSE,
  fish = FALSE
)

myResDN_fe

# compare this to the other results
myResDN_degree_2

myResDN_cont

# as you can see, degree controls make quite a difference IN THIS CASE,
# but which one is used is less important

# comparison to without degree controls
myEffects_cont <- createEffectsObject(
  list(
    list("loops"),
    list("reciprocity_min"),
    list("dyadic_covariate", attribute.index = "sameRegion")
  )
)

myAlg_cont <- createAlgorithm(myState_fe, myEffects_cont, multinomialProposal = FALSE,
                            initGain = 0.8)

myResDN_cont <- estimateMobilityNetwork(
  myState_fe, myCache_fe, myEffects_cont, myAlg_cont,
  parallel = TRUE, cpus = 4,
  verbose = TRUE,
  returnDeps = FALSE,
  fish = FALSE
)

myResDN_cont






