estimate_logL {ERPM}

Estimate log likelihood

Description

Function to estimate the log likelihood of a model for an observed partition

Usage

estimate_logL(
  partition,
  nodes,
  effects,
  objects,
  theta,
  theta_0,
  M,
  num.steps,
  burnin,
  thining,
  neighborhoods = c(0.7, 0.3, 0),
  numgroups.allowed = NULL,
  numgroups.simulated = NULL,
  sizes.allowed = NULL,
  sizes.simulated = NULL,
  logL_0 = NULL,
  parallel = FALSE,
  cpus = 1,
  verbose = FALSE
)

Arguments

partition

observed partition

nodes

node set (data frame)

effects

effects/sufficient statistics (list with a vector "names", and a vector "objects")

objects

objects used for statistics calculation (list with a vector "name", and a vector "object")

theta

estimated model parameters

theta_0

model parameters if all other effects than "num-groups" are fixed to 0 (basic Dirichlet partition model)

M

number of steps in the path-sampling algorithm

num.steps

number of samples in each step

burnin

integer for the number of burn-in steps before sampling

thining

integer for the number of thining steps between sampling

neighborhoods

= c(0.7,0.3,0) way of choosing partitions

numgroups.allowed

= NULL, # vector containing the number of groups allowed in the partition (now, it only works with vectors like num_min:num_max)

numgroups.simulated

= NULL, # vector containing the number of groups simulated

sizes.allowed

= NULL, vector of group sizes allowed in sampling (now, it only works for vectors like size_min:size_max)

sizes.simulated

= NULL, vector of group sizes allowed in the Markov chain but not necessraily sampled (now, it only works for vectors like size_min:size_max)

logL_0

= NULL, if known, the value of the log likelihood of the basic dirichlet model

parallel

= FALSE, indicating whether the code should be run in parallel

cpus

= 1, number of cpus required for the parallelization

verbose

= FALSE, to print the current step the algorithm is in

Value

List with the log likelihood , AIC, lambda and the draws

Examples

# estimate the log-likelihood and AIC of an estimated model (e.g. useful to compare two models)

# define an arbitrary set of n = 6 nodes with attributes, and an arbitrary covariate matrix
n <- 6
nodes <- data.frame(label = c("A","B","C","D","E","F"),
                    gender = c(1,1,2,1,2,2),
                    age = c(20,22,25,30,30,31)) 
friendship <- matrix(c(0, 1, 1, 1, 0, 0,
                       1, 0, 0, 0, 1, 0,
                       1, 0, 0, 0, 1, 0,
                       1, 0, 0, 0, 0, 0,
                       0, 1, 1, 0, 0, 1,
                       0, 0, 0, 0, 1, 0), 6, 6, TRUE)

# choose the effects to be included (see manual for all effect names)
effects <- list(names = c("num_groups","same","diff","tie"),
objects = c("partition","gender","age","friendship"))
objects <- list()
objects[[1]] <- list(name = "friendship", object = friendship)

# define observed partition 
partition <- c(1,1,2,2,2,3)
# (an exemplary estimation is internally stored in order to save time)

# first: estimate the ML estimates of a simple model with only one parameter 
# for number of groups (this parameter should be in the model!)
likelihood_function <- function(x){ exp(x*max(partition)) / compute_numgroups_denominator(n,x)}
curve(likelihood_function, from=-2, to=0)
parameter_base <- optimize(likelihood_function, interval=c(-2, 0), maximum=TRUE)
parameters_basemodel <- c(parameter_base$maximum,0,0,0)


# estimate logL and AIC
logL_AIC <- estimate_logL(partition,
                          nodes,
                          effects, 
                          objects,
                          theta = estimation$results$est,
                          theta_0 = parameters_basemodel,
                          M = 3,
                          num.steps = 200,
                          burnin = 100,
                          thining = 20)
logL_AIC$logL
logL_AIC$AIC