Reference

ALRsimple

An autologistic regression model with "simple smoothing": the unary parameter is of type LinPredUnary, and the pairwise parameter is of type SimplePairwise.

Constructors

ALRsimple(unary::LinPredUnary, pairwise::SimplePairwise;
    Y::Union{Nothing,<:VecOrMat} = nothing,
    centering::CenteringKinds = none, 
    coding::Tuple{Real,Real} = (-1,1),
    labels::Tuple{String,String} = ("low","high"), 
    coordinates::SpatialCoordinates = [(0.0,0.0) for i=1:size(unary,1)]
)
ALRsimple(graph::SimpleGraph{Int}, X::Float2D3D; 
    Y::VecOrMat = Array{Bool,2}(undef,nv(graph),size(X,3)),
    β::Vector{Float64} = zeros(size(X,2)),
    λ::Float64 = 0.0, 
    centering::CenteringKinds = none, 
    coding::Tuple{Real,Real} = (-1,1),
    labels::Tuple{String,String} = ("low","high"),
    coordinates::SpatialCoordinates = [(0.0,0.0) for i=1:nv(graph)]
)

Arguments

• Y: the array of dichotomous responses. Any array with 2 unique values will work. If the array has only one unique value, it must equal one of the coding values. The supplied object will be internally represented as a Boolean array.
• β: the regression coefficients.
• λ: the association parameter.
• centering: controls what form of centering to use.
• coding: determines the numeric coding of the dichotomous responses.
• labels: a 2-tuple of text labels describing the meaning of Y. The first element is the label corresponding to the lower coding value.
• coordinates: an array of 2- or 3-tuples giving spatial coordinates of each vertex in the graph.

Examples

julia> using Graphs
julia> X = rand(10,3);            #-predictors
julia> Y = rand([-2, 3], 10);     #-responses
julia> g = Graph(10,20);          #-graph
julia> u = LinPredUnary(X);
julia> p = SimplePairwise(g);
julia> model1 = ALRsimple(u, p, Y=Y);
julia> model2 = ALRsimple(g, X, Y=Y);
julia> all([getfield(model1, fn)==getfield(model2, fn) for fn in fieldnames(ALRsimple)])
true

ALfit

A type to hold estimation output for autologistic models. Fitting functions return an object of this type.

Depending on the fitting method, some fields might not be set. Fields that are not used are set to nothing or to zero-dimensional arrays. The fields are:

• estimate: A vector of parameter estimates.
• se: A vector of standard errors for the estimates.
• pvalues: A vector of p-values for testing the null hypothesis that the parameters equal zero (one-at-a time hypothesis tests).
• CIs: A vector of 95% confidence intervals for the parameters (a vector of 2-tuples).
• optim: the output of the call to optimize used to get the estimates.
• Hinv (used by fit_ml!): The inverse of the Hessian matrix of the objective function, evaluated at the estimate.
• nboot (fit_pl!): number of bootstrap samples to use for error estimation.
• kwargs (fit_pl!): holds extra keyword arguments passed in the call to the fitting function.
• bootsamples (fit_pl!): the bootstrap samples.
• bootestimates (fit_pl!): the bootstrap parameter estimates.
• convergence: either a Boolean indicating optimization convergence ( for fit_ml!), or a vector of such values for the optimizations done to estimate bootstrap replicates.

The empty constructor ALfit() will initialize an object with all fields empty, so the needed fields can be filled afterwards.

ALfull

An autologistic model with a FullUnary unary parameter type and a FullPairwise pairwise parameter type. This model has the maximum number of unary parameters (one parameter per variable per observation), and an association matrix with one parameter per edge in the graph.

Constructors

ALfull(unary::FullUnary, pairwise::FullPairwise; 
    Y::Union{Nothing,<:VecOrMat}=nothing, 
    centering::CenteringKinds=none, 
    coding::Tuple{Real,Real}=(-1,1),
    labels::Tuple{String,String}=("low","high"), 
    coordinates::SpatialCoordinates=[(0.0,0.0) for i=1:size(unary,1)]
)
ALfull(graph::SimpleGraph{Int}, alpha::Float1D2D, lambda::Vector{Float64}; 
    Y::VecOrMat=Array{Bool,2}(undef,nv(graph),size(alpha,2)), 
    centering::CenteringKinds=none, 
    coding::Tuple{Real,Real}=(-1,1),
    labels::Tuple{String,String}=("low","high"),
    coordinates::SpatialCoordinates=[(0.0,0.0) for i=1:nv(graph)]
)
ALfull(graph::SimpleGraph{Int}, count::Int=1; 
    Y::VecOrMat=Array{Bool,2}(undef,nv(graph),count), 
    centering::CenteringKinds=none, 
    coding::Tuple{Real,Real}=(-1,1),
    labels::Tuple{String,String}=("low","high"),
    coordinates::SpatialCoordinates=[(0.0,0.0) for i=1:nv(graph)]
)

Arguments

• Y: the array of dichotomous responses. Any array with 2 unique values will work. If the array has only one unique value, it must equal one of th coding values. The supplied object will be internally represented as a Boolean array.
• centering: controls what form of centering to use.
• coding: determines the numeric coding of the dichotomous responses.
• labels: a 2-tuple of text labels describing the meaning of Y. The first element is the label corresponding to the lower coding value.
• coordinates: an array of 2- or 3-tuples giving spatial coordinates of each vertex in the graph.

Examples

julia> g = Graph(10, 20);          #-graph (20 edges)
julia> alpha = zeros(10, 4);       #-unary parameter values
julia> lambda = rand(20);          #-pairwise parameter values
julia> Y = rand([0, 1], 10, 4);    #-responses
julia> u = FullUnary(alpha);
julia> p = FullPairwise(g, 4);
julia> setparameters!(p, lambda);
julia> model1 = ALfull(u, p, Y=Y);
julia> model2 = ALfull(g, alpha, lambda, Y=Y);
julia> model3 = ALfull(g, 4, Y=Y);
julia> setparameters!(model3, [alpha[:]; lambda]);
julia> all([getfield(model1, fn)==getfield(model2, fn)==getfield(model3, fn)
            for fn in fieldnames(ALfull)])
true

ALsimple

An autologistic model with a FullUnary unary parameter type and a SimplePairwise pairwise parameter type. This model has the maximum number of unary parameters (one parameter per variable per observation), and a single association parameter.

Constructors

ALsimple(unary::FullUnary, pairwise::SimplePairwise; 
    Y::Union{Nothing,<:VecOrMat} = nothing, 
    centering::CenteringKinds = none, 
    coding::Tuple{Real,Real} = (-1,1),
    labels::Tuple{String,String} = ("low","high"), 
    coordinates::SpatialCoordinates = [(0.0,0.0) for i=1:size(unary,1)]
)
ALsimple(graph::SimpleGraph{Int}, alpha::Float1D2D; 
    Y::VecOrMat = Array{Bool,2}(undef,nv(graph),size(alpha,2)), 
    λ::Float64 = 0.0, 
    centering::CenteringKinds = none, 
    coding::Tuple{Real,Real} = (-1,1),
    labels::Tuple{String,String} = ("low","high"),
    coordinates::SpatialCoordinates = [(0.0,0.0) for i=1:nv(graph)]
)
ALsimple(graph::SimpleGraph{Int}, count::Int = 1; 
    Y::VecOrMat = Array{Bool,2}(undef,nv(graph),size(alpha,2)), 
    λ::Float64=0.0, 
    centering::CenteringKinds=none, 
    coding::Tuple{Real,Real}=(-1,1),
    labels::Tuple{String,String}=("low","high"),
    coordinates::SpatialCoordinates=[(0.0,0.0) for i=1:nv(graph)]
)

Arguments

• Y: the array of dichotomous responses. Any array with 2 unique values will work. If the array has only one unique value, it must equal one of th coding values. The supplied object will be internally represented as a Boolean array.
• λ: the association parameter.
• centering: controls what form of centering to use.
• coding: determines the numeric coding of the dichotomous responses.
• labels: a 2-tuple of text labels describing the meaning of Y. The first element is the label corresponding to the lower coding value.
• coordinates: an array of 2- or 3-tuples giving spatial coordinates of each vertex in the graph.

Examples

julia> alpha = zeros(10, 4);       #-unary parameter values
julia> Y = rand([0, 1], 10, 4);    #-responses
julia> g = Graph(10, 20);          #-graph
julia> u = FullUnary(alpha);
julia> p = SimplePairwise(g, 4);
julia> model1 = ALsimple(u, p, Y=Y);
julia> model2 = ALsimple(g, alpha, Y=Y);
julia> model3 = ALsimple(g, 4, Y=Y);
julia> all([getfield(model1, fn)==getfield(model2, fn)==getfield(model3, fn)
            for fn in fieldnames(ALsimple)])
true

AbstractAutologisticModel

Abstract type representing autologistic models. All concrete subtypes should have the following fields:

• responses::Array{Bool,2} – The binary observations. Rows are for nodes in the graph, and columns are for independent (vector) observations. It is a 2D array even if there is only one observation.
• unary<:AbstractUnaryParameter – Specifies the unary part of the model.
• pairwise<:AbstractPairwiseParameter – Specifies the pairwise part of the model (including the graph).
• centering<:CenteringKinds – Specifies the form of centering used, if any.
• coding::Tuple{T,T} where T<:Real – Gives the numeric coding of the responses.
• labels::Tuple{String,String} – Provides names for the high and low states.
• coordinates<:SpatialCoordinates – Provides 2D or 3D coordinates for each vertex in the graph.

This type has the following functions defined, considered part of the type's interface. They cover most operations one will want to perform. Concrete subtypes should not have to define custom overrides unless more specialized or efficient algorithms exist for the subtype.

• getparameters and setparameters!
• getunaryparameters and setunaryparameters!
• getpairwiseparameters and setpairwiseparameters!
• centeringterms
• negpotential, pseudolikelihood, and loglikelihood
• fullPMF, marginalprobabilities, and conditionalprobabilities
• fit_pl! and fit_ml!
• sample and oneboot
• showfields

Examples

julia> M = ALsimple(Graph(4,4));
julia> typeof(M)
ALsimple{CenteringKinds,Int64,Nothing}
julia> isa(M, AbstractAutologisticModel)
true

AbstractPairwiseParameter

Abstract type representing the pairwise part of an autologistic regression model.

All concrete subtypes should have the following fields:

• G::SimpleGraph{Int} – The graph for the model.
• count::Int – The number of observations.

In addition to getindex() and setindex!(), any concrete subtype P<:AbstractPairwiseParameter should also have the following methods defined:

• getparameters(P), returning a Vector{Float64}
• setparameters!(P, newpar::Vector{Float64}) for setting parameter values.

Note that indexing is performance-critical and should be implemented carefully in subtypes.

The intention is that each subtype should implement a different way of parameterizing the association matrix. The way parameters are stored and values computed is up to the subtypes.

This type inherits from AbstractArray{Float64, 3}. The third index is to allow for multiple observations. P[:,:,r] should return the association matrix of the rth observation in an appropriate subtype of AbstractMatrix. It is not intended that the third index will be used for range or vector indexing like P[:,:,1:5] (though this may work due to AbstractArray fallbacks).

Examples

julia> M = ALsimple(Graph(4,4));
julia> typeof(M.pairwise)
SimplePairwise
julia> isa(M.pairwise, AbstractPairwiseParameter)
true

AbstractUnaryParameter

Abstract type representing the unary part of an autologistic regression model.

This type inherits from AbstractArray{Float64, 2}. The first dimension is for vertices/variables in the graph, and the second dimension is for observations. It is two-dimensional even if there is only one observation.

Implementation details are left to concrete subtypes, and will depend on how the unary terms are parametrized. Note that indexing is performance-critical.

Concrete subtypes should implement getparameters, setparameters!, and showfields.

Examples

julia> M = ALsimple(Graph(4,4));
julia> typeof(M.unary)
FullUnary
julia> isa(M.unary, AbstractUnaryParameter)
true

CenteringKinds

An enumeration to facilitate choosing a form of centering for the model. Available choices are:

• none: no centering (centering adjustment equals zero).
• expectation: the centering adjustment is the expected value of the response under the assumption of independence (this is what has been used in the "centered autologistic model").
• onehalf: a constant value of centering adjustment equal to 0.5 (this produces the "symmetric autologistic model" when used with 0,1 coding).

The default/recommended model has centering of none with (-1, 1) coding.

Examples

julia> CenteringKinds
Enum CenteringKinds:
none = 0
expectation = 1
onehalf = 2

FullPairwise

A type representing an association matrix with a "saturated" parametrization–one parameter for each edge in the graph.

In this type, the association matrix for each observation is a symmetric matrix with the same pattern of nonzeros as the graph's adjacency matrix, but with arbitrary values in those locations. The package convention is to provide parameters as a vector of Float64. So getparameters and setparameters! use a vector of ne(G) values that correspond to the nonzero locations in the upper triangle of the adjacency matrix, in the same (lexicographic) order as edges(G).

The association matrix is stored as a SparseMatrixCSC{Float64,Int64} in the field Λ.

This type does not allow for different observations to have different association matricse. So while size returns a 3-dimensional result, the third index is ignored when accessing the array's elements.

Constructors

FullPairwise(G::SimpleGraph, count::Int=1)
FullPairwise(n::Int, count::Int=1)
FullPairwise(λ::Real, G::SimpleGraph)
FullPairwise(λ::Real, G::SimpleGraph, count::Int)
FullPairwise(λ::Vector{Float64}, G::SimpleGraph)

If provide only a graph, set λ = zeros(nv(G)). If provide only an integer, set λ to zeros and make a totally disconnected graph. If provide a graph and a scalar, convert the scalar to a vector of the right length.

Examples

julia> g = makegrid4(2,2).G;
julia> λ = [1.0, 2.0, -1.0, -2.0];
julia> p = FullPairwise(λ, g);
julia> typeof(p.Λ)
SparseArrays.SparseMatrixCSC{Float64,Int64}

julia> Matrix(p[:,:])
4×4 Array{Float64,2}:
 0.0   1.0   2.0   0.0
 1.0   0.0   0.0  -1.0
 2.0   0.0   0.0  -2.0
 0.0  -1.0  -2.0   0.0

FullUnary

The unary part of an autologistic model, with one parameter per vertex per observation. The type has only a single field, for holding an array of parameters.

Constructors

FullUnary(alpha::Array{Float64,1}) 
FullUnary(n::Int)                     #-initializes parameters to zeros
FullUnary(n::Int, m::Int)             #-initializes parameters to zeros

Examples

julia> u = FullUnary(5, 3);
julia> u[:,:]
5×3 Array{Float64,2}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

LinPredUnary

The unary part of an autologistic model, parametrized as a regression linear predictor. Its fields are X, an n-by-p-by-m matrix (n obs, p predictors, m observations), and β, a p-vector of parameters (the same for all observations).

Constructors

LinPredUnary(X::Matrix{Float64}, β::Vector{Float64})
LinPredUnary(X::Matrix{Float64})
LinPredUnary(X::Array{Float64, 3})
LinPredUnary(n::Int,p::Int)
LinPredUnary(n::Int,p::Int,m::Int)

Any quantities not provided in the constructors are initialized to zeros.

Examples

julia> u = LinPredUnary(ones(5,3,2), [1.0, 2.0, 3.0]);
julia> u[:,:]
5×2 Array{Float64,2}:
 6.0  6.0
 6.0  6.0
 6.0  6.0
 6.0  6.0
 6.0  6.0

SamplingMethods

An enumeration to facilitate choosing a method for random sampling from autologistic models. Available choices are:

• Gibbs: Gibbs sampling.
• perfect_bounding_chain: Perfect sampling, using a bounding chain algorithm.
• perfect_reuse_samples: Perfect sampling. CFTP implemented by reusing random numbers.
• perfect_reuse_seeds: Perfect sampling. CFTP implemented by reusing RNG seeds.
• perfect_read_once: Perfect sampling. Read-once CFTP implementation.

All of the perfect sampling methods are implementations of coupling from the past (CFTP). perfect_bounding_chain uses a bounding chain approach that holds even when Λ contains negative elements; the other three options rely on a monotonicity argument that requires Λ to have only positive elements (though they should work similar to Gibbs sampling in that case).

Different perfect sampling implementations might work best for different models, and parameter settings exist where perfect sampling coalescence might take a prohibitively long time. For these reasons, Gibbs sampling is the default in sample.

Examples

julia> SamplingMethods
Enum SamplingMethods:
Gibbs = 0
perfect_reuse_samples = 1
perfect_reuse_seeds = 2
perfect_read_once = 3
perfect_bounding_chain = 4

SimplePairwise

Pairwise association matrix, parametrized as a scalar parameter times the adjacency matrix.

Constructors

SimplePairwise(G::SimpleGraph, count::Int=1) SimplePairwise(n::Int, count::Int=1) SimplePairwise(λ::Real, G::SimpleGraph) SimplePairwise(λ::Real, G::SimpleGraph, count::Int)

If provide only a graph, set λ = 0. If provide only an integer, set λ = 0 and make a totally disconnected graph. If provide a graph and a scalar, convert the scalar to a length-1 vector.

Every observation must have the same association matrix in this case. So while we internally treat it like an n-by-n-by-m matrix, just return a 2D n-by-n matrix to the user.

Examples

julia> g = makegrid4(2,2).G;
julia> λ = 1.0;
julia> p = SimplePairwise(λ, g, 4);    #-4 observations
julia> size(p)
(4, 4, 4)

julia> Matrix(p[:,:,:])
4×4 Array{Float64,2}:
 0.0  1.0  1.0  0.0
 1.0  0.0  0.0  1.0
 1.0  0.0  0.0  1.0
 0.0  1.0  1.0  0.0

Type alias for Union{Array{Float64,1},Array{Float64,2}}

Type alias for Union{Array{Float64,2},Array{Float64,3}}

Type alias for Union{Array{NTuple{2,T},1},Array{NTuple{3,T},1}} where T<:Real

Type alias for Union{Array{T,1}, Array{T,2}} where T

addboot!(fit::ALfit, bootsamples::Array{Float64,3}, 
         bootestimates::Array{Float64,2}, convergence::Vector{Bool})

Add parametric bootstrap information in arrays bootsamples, bootestimates, and convergence to model fitting information fit. If fit already contains bootstrap data, the new data is appended to the existing data, and statistics are recomputed.

Examples

julia> using Random;
julia> Random.seed!(1234);
julia> G = makegrid4(4,3).G;
julia> Y=[[fill(-1,4); fill(1,8)] [fill(-1,3); fill(1,9)] [fill(-1,5); fill(1,7)]];
julia> model = ALRsimple(G, ones(12,1,3), Y=Y);
julia> fit = fit_pl!(model, start=[-0.4, 1.1]);
julia> samps = zeros(12,3,10);
julia> ests = zeros(2,10);
julia> convs = fill(false, 10);
julia> for i = 1:10
           temp = oneboot(model, start=[-0.4, 1.1])
           samps[:,:,i] = temp.sample
           ests[:,i] = temp.estimate
           convs[i] = temp.convergence
       end
julia> addboot!(fit, samps, ests, convs)
julia> summary(fit)
name          est     se      95% CI
parameter 1   -0.39   0.442      (-1.09, 0.263)
parameter 2    1.1    0.279   (-0.00664, 0.84)

addboot!(fit::ALfit, bootresults::Array{T,1}) where 
    T <: NamedTuple{(:sample, :estimate, :convergence)}

An addboot! method taking bootstrap data as an array of named tuples. Tuples are of the form produced by oneboot.

Examples

julia>     using Random;
julia>     Random.seed!(1234);
julia>     G = makegrid4(4,3).G;
julia>     Y=[[fill(-1,4); fill(1,8)] [fill(-1,3); fill(1,9)] [fill(-1,5); fill(1,7)]];
julia>     model = ALRsimple(G, ones(12,1,3), Y=Y);
julia>     fit = fit_pl!(model, start=[-0.4, 1.1]);
julia>     boots = [oneboot(model, start=[-0.4, 1.1]) for i = 1:10];
julia>     addboot!(fit, boots)
julia>     summary(fit)
name          est     se      95% CI
parameter 1   -0.39   0.442      (-1.09, 0.263)
parameter 2    1.1    0.279   (-0.00664, 0.84)

centeringterms(M::AbstractAutologisticModel, kind::Union{Nothing,CenteringKinds}=nothing)

Returns an Array{Float64,2} of the same dimension as M.unary, giving the centering adjustments for autologistic model M. centeringterms(M,kind) returns the centering adjustment that would be used if centering were of type kind.

Examples

julia> G = makegrid8(2,2).G;
julia> X = [ones(4) [-2; -1; 1; 2]];
julia> M1 = ALRsimple(G, X, β=[-1.0, 2.0]);                 #-No centering (default)
julia> M2 = ALRsimple(G, X, β=[-1.0, 2.0], centering=expectation);  #-Centered model
julia> [centeringterms(M1) centeringterms(M2) centeringterms(M1, onehalf)]
4×3 Array{Float64,2}:
 0.0  -0.999909  0.5
 0.0  -0.995055  0.5
 0.0   0.761594  0.5
 0.0   0.995055  0.5

conditionalprobabilities(M::AbstractAutologisticModel; vertices=1:size(M.unary)[1], 
                         indices=1:size(M.unary,2))

Compute the conditional probability that variables in autologistic model M take the high state, given the current values of all of their neighbors. If vertices or indices are provided, the results are only computed for the desired variables & observations. Otherwise results are computed for all variables and observations.

Examples

julia> Y = [ones(9) zeros(9)];
julia> G = makegrid4(3,3).G;
julia> model = ALsimple(G, ones(9,2), Y=Y, λ=0.5);    #-Variables on a 3×3 grid, 2 obs.
julia> conditionalprobabilities(model, vertices=5)    #-Cond. probs. of center vertex.
1×2 Array{Float64,2}:
 0.997527  0.119203

julia> conditionalprobabilities(model, indices=2)     #-Cond probs, 2nd observation.
9×1 Array{Float64,2}:
 0.5
 0.26894142136999516
 0.5
 0.26894142136999516
 0.11920292202211756
 0.26894142136999516
 0.5
 0.26894142136999516
 0.5

`Autologistic.datasets(name::String)`

Open data sets for demonstrations of Autologistic regression. Returns the data set as a DataFrame.

fit_ml!(M::AbstractAutologisticModel;
    start=zeros(length(getparameters(M))),
    verbose::Bool=false,
    force::Bool=false,
    kwargs...
)

Fit autologistic model M using maximum likelihood. Will fail for models with more than 20 vertices, unless force=true. Use fit_pl! for larger models.

Arguments

• start: initial value to use for optimization.
• verbose: should progress information be printed to the console?
• force: set to true to force computation of the likelihood for large models.
• kwargs... extra keyword arguments that are passed on to optimize().

Examples

julia> G = makegrid4(4,3).G;
julia> model = ALRsimple(G, ones(12,1), Y=[fill(-1,4); fill(1,8)]);
julia> mle = fit_ml!(model);
julia> summary(mle)
name          est      se      p-value   95% CI
parameter 1   0.0791   0.163   0.628       (-0.241, 0.399)
parameter 2   0.425    0.218   0.0511    (-0.00208, 0.852)

fit_pl!(M::AbstractAutologisticModel;
    start=zeros(length(getparameters(M))), 
    verbose::Bool=false,
    nboot::Int = 0,
    kwargs...)

Fit autologistic model M using maximum pseudolikelihood.

Arguments

• start: initial value to use for optimization.
• verbose: should progress information be printed to the console?
• nboot: number of samples to use for parametric bootstrap error estimation. If nboot=0 (the default), no bootstrap is run.
• kwargs... extra keyword arguments that are passed on to optimize() or sample(), as appropriate.

Examples

julia> Y=[[fill(-1,4); fill(1,8)] [fill(-1,3); fill(1,9)] [fill(-1,5); fill(1,7)]];
julia> model = ALRsimple(makegrid4(4,3).G, ones(12,1,3), Y=Y);
julia> fit = fit_pl!(model, start=[-0.4, 1.1]);
julia> summary(fit)
name          est     se   p-value   95% CI
parameter 1   -0.39
parameter 2    1.1

fullPMF(M::AbstractAutologisticModel; 
    indices=1:size(M.unary,2), 
    force::Bool=false
)

Compute the PMF of an AbstractAutologisticModel, and return a NamedTuple (:table, :partition).

For an AbstractAutologisticModel with $n$ variables and $m$ observations, :table is a $2^n×(n+1)×m$ array of Float64. Each page of the 3D array holds a probability table for an observation. Each row of the table holds a specific configuration of the responses, with the corresponding probability in the last column. In the $m=1$ case, :table is a 2D array.

Output :partition is a vector of normalizing constant (a.k.a. partition function) values. In the $m=1$ case, it is a scalar Float64.

Arguments

• indices: indices of specific observations from which to obtain the output. By default, all observations are used.
• force: calling the function with $n>20$ will throw an error unless force=true.

Examples

julia> M = ALRsimple(Graph(3,0),ones(3,1));
julia> pmf = fullPMF(M);
julia> pmf.table
8×4 Array{Float64,2}:
 -1.0  -1.0  -1.0  0.125
 -1.0  -1.0   1.0  0.125
 -1.0   1.0  -1.0  0.125
 -1.0   1.0   1.0  0.125
  1.0  -1.0  -1.0  0.125
  1.0  -1.0   1.0  0.125
  1.0   1.0  -1.0  0.125
  1.0   1.0   1.0  0.125
julia> pmf.partition
 8.0

getpairwiseparameters(M::AbstractAutologisticModel)

Extracts the pairwise parameters from an autologistic model. Parameters are always returned as an Array{Float64,1}.

getparameters(x)

A generic function for extracting the parameters from an autologistic model, a unary term, or a pairwise term. Parameters are always returned as an Array{Float64,1}. If typeof(x) <: AbstractAutologisticModel, the returned vector is partitioned with the unary parameters first.

getunaryparameters(M::AbstractAutologisticModel)

Extracts the unary parameters from an autologistic model. Parameters are always returned as an Array{Float64,1}.

loglikelihood(M::AbstractAutologisticModel; 
    force::Bool=false
)

Compute the natural logarithm of the likelihood for autologistic model M. This will throw an error for models with more than 20 vertices, unless force=true.

Examples

julia> model = ALRsimple(makegrid4(2,2)[1], ones(4,2,3), centering = expectation,
                         coding = (0,1), Y = repeat([true, true, false, false],1,3));
julia> setparameters!(model, [1.0, 1.0, 1.0]);
julia> loglikelihood(model)
-11.86986109487605

makebool(v::VecOrMat, vals=nothing)

Makes a 2D array of Booleans out of a 1- or 2-D input. The 2nd argument vals optionally can be a 2-tuple (low, high) specifying the two possible values in v (useful for the case where all elements of v take one value or the other).

• If v has more than 2 unique values, throws an error.
• If v has exactly 2 unique values, use those to set the coding (ignore vals).
• If v has 1 unique value, use vals to determine if it's the high or low value (throw an error if the single value isn't in vals).

Examples

julia> makebool([1.0 2.0; 1.0 2.0])
2×2 Array{Bool,2}:
 false  true
 false  true

julia> makebool(["yes", "no", "no"])
3×1 Array{Bool,2}:
  true
 false
 false

julia> [makebool([1, 1, 1], (-1,1)) makebool([1, 1, 1], (1, 2))]
3×2 Array{Bool,2}:
 true  false
 true  false
 true  false

makecoded(M::AbstractAutologisticModel)

A convenience method for makecoded(M.responses, M.coding). Use it to retrieve a model's responses in coded form.

Examples

julia> M = ALRsimple(Graph(4,3), rand(4,2), Y=[true, false, false, true], coding=(-1,1));
julia> makecoded(M)
4×1 Array{Float64,2}:
  1.0
 -1.0
 -1.0
  1.0

makecoded(b::VecOrMat, coding::Tuple{Real,Real})

Convert Boolean responses into coded values. The first argument is boolean. Returns a 2D array of Float64.

Examples

julia> makecoded([true, false, false, true], (-1, 1))
4×1 Array{Float64,2}:
  1.0
 -1.0
 -1.0
  1.0

makegrid4(r::Int, c::Int, xlim::Tuple{Real,Real}=(0.0,1.0), 
          ylim::Tuple{Real,Real}=(0.0,1.0))

Returns a named tuple (:G, :locs), where :G is a graph, and :locs is an array of numeric tuples. Vertices of :G are laid out in a rectangular, 4-connected grid with r rows and c columns. The tuples in :locs contain the spatial coordinates of each vertex. Optional arguments xlim and ylim determine the bounds of the rectangular layout.

Examples

julia> out4 = makegrid4(11, 21, (-1,1), (-10,10));
julia> nv(out4.G) == 11*21                  #231
true
julia> ne(out4.G) == 11*20 + 21*10          #430
true
julia> out4.locs[11*10 + 6] == (0.0, 0.0)   #location of center vertex.
true

makegrid8(r::Int, c::Int, xlim::Tuple{Real,Real}=(0.0,1.0), 
          ylim::Tuple{Real,Real}=(0.0,1.0))

Returns a named tuple (:G, :locs), where :G is a graph, and :locs is an array of numeric tuples. Vertices of :G are laid out in a rectangular, 8-connected grid with r rows and c columns. The tuples in :locs contain the spatial coordinates of each vertex. Optional arguments xlim and ylim determine the bounds of the rectangular layout.

Examples

julia> out8 = makegrid8(11, 21, (-1,1), (-10,10));
julia> nv(out8.G) == 11*21                      #231
true
julia> ne(out8.G) == 11*20 + 21*10 + 2*20*10    #830
true
julia> out8.locs[11*10 + 6] == (0.0, 0.0)       #location of center vertex.
true

makespatialgraph(coords::C, δ::Real) where C<:SpatialCoordinates

Returns a named tuple (:G, :locs), where :G is a graph, and :locs is an array of numeric tuples. Each element of coords is a 2- or 3-tuple of spatial coordinates, and this argument is returned unchanged as :locs. The graph :G has length(coords) vertices, with edges connecting every pair of vertices within Euclidean distance δ of each other.

Examples

julia> c = [(Float64(i), Float64(j)) for i = 1:5 for j = 1:5];
julia> out = makespatialgraph(c, sqrt(2));
julia> out.G
{25, 72} undirected simple Int64 graph

julia> length(out.locs)
25

marginalprobabilities(M::AbstractAutologisticModel;
    indices=1:size(M.unary,2), 
    force::Bool=false
)

Compute the marginal probability that variables in autologistic model M takes the high state. For a model with n vertices and m observations, returns an n-by-m array (or an n-vector if m==1). The [i,j]th element is the marginal probability of the high state in the ith variable at the jth observation.

This function computes the exact marginals. For large models, approximate the marginal probabilities by sampling, e.g. sample(M, ..., average=true).

Arguments

• indices: used to return only the probabilities for certain observations.
• force: the function will throw an error for n > 20 unless force=true.

Examples

julia> M = ALsimple(Graph(3,0), [[-1.0; 0.0; 1.0] [-1.0; 0.0; 1.0]])
julia> marginalprobabilities(M)
3×2 Array{Float64,2}:
 0.119203  0.119203
 0.5       0.5
 0.880797  0.880797

negpotential(M::AbstractAutologisticModel)

Returns an m-vector of Float64 negpotential values, where m is the number of observations found in M.responses.

Examples

julia> M = ALsimple(makegrid4(3,3).G, ones(9));
julia> f = fullPMF(M);
julia> exp(negpotential(M)[1])/f.partition ≈ exp(loglikelihood(M))
true

oneboot(M::AbstractAutologisticModel, params::Vector{Float64};
    start=zeros(length(getparameters(M))),
    verbose::Bool=false,
    kwargs...
)

Computes one bootstrap replicate using model M, but using parameters params for generating samples, instead of getparameters(M).

oneboot(M::AbstractAutologisticModel; 
    start=zeros(length(getparameters(M))),
    verbose::Bool=false,
    kwargs...
)

Performs one parametric bootstrap replication from autologistic model M: draw a random sample from M, use that sample as the responses, and re-fit the model. Returns a named tuple (:sample, :estimate, :convergence), where :sample holds the random sample, :estimate holds the parameter estimates, and :convergence holds a bool indicating whether or not the optimization converged. The parameters of M remain unchanged by calling oneboot.

Arguments

• start: starting parameter values to use for optimization
• verbose: should progress information be written to the console?
• kwargs...: extra keyword arguments that are passed to optimize() or sample(), as appropriate.

Examples

jldoctest julia> G = makegrid4(4,3).G; julia> model = ALRsimple(G, ones(12,1), Y=[fill(-1,4); fill(1,8)]); julia> theboot = oneboot(model, method=Gibbs, burnin=250); julia> fieldnames(typeof(theboot)) (:sample, :estimate, :convergence)

pseudolikelihood(M::AbstractAutologisticModel)

Computes the negative log pseudolikelihood for autologistic model M. Returns a Float64.

Examples

julia> X = [1.1 2.2
            1.0 2.0
            2.1 1.2
            3.0 0.3];
julia> Y = [0; 0; 1; 0];
julia> M3 = ALRsimple(makegrid4(2,2)[1], cat(X,X,dims=3), Y=cat(Y,Y,dims=2), 
                      β=[-0.5, 1.5], λ=1.25, centering=expectation);
julia> pseudolikelihood(M3)
12.333549445795818

setpairwiseparameters!(M::AbstractAutologisticModel, newpars::Vector{Float64})

Sets the pairwise parameters of autologistic model M to the values in newpars.

setparameters!(x, newpars::Vector{Float64})

A generic function for setting the parameter values of an autologistic model, a unary term, or a pairwise term. Parameters are always passed as an Array{Float64,1}. If typeof(x) <: AbstractAutologisticModel, the newpars is assumed partitioned with the unary parameters first.

setunaryparameters!(M::AbstractAutologisticModel, newpars::Vector{Float64})

Sets the unary parameters of autologistic model M to the values in newpars.

sample(M::AbstractAutologisticModel, k::Int = 1;
    method::SamplingMethods = Gibbs,
    indices = 1:size(M.unary,2), 
    average::Bool = false, 
    config = nothing, 
    burnin::Int = 0,
    skip::Int = 0,
    verbose::Bool = false
)

Draws k random samples from autologistic model M. For a model with n vertices in its graph, the return value is:

• When average=false, an n × length(indices) × k array, with singleton dimensions dropped. This array holds the random samples.
• When average=true, an n × length(indices) array, with singleton dimensions dropped. This array holds the estimated marginal probabilities of observing the "high" level at each vertex.

Arguments

• method: a member of the enum SamplingMethods, specifying which sampling method will be used. The default is Gibbs sampling. Where feasible, it is recommended to use one of the perfect sampling alternatives. See SamplingMethods for more.
• indices: gives the indices of the observation to use for sampling. If the model has more than one observation, then k samples are drawn for each observation's parameter settings. Use indices to restrict the samples to a subset of observations.
• average: controls the form of the output. When average=true, the return value is the proportion of "high" samples at each vertex. (Note that this is not actually the arithmetic average of the samples, unless the coding is (0,1). Rather, it is an estimate of the probability of getting a "high" outcome.) When average=false, the full set of samples is returned.
• config: allows a starting configuration of the random variables to be provided. Only used if method=Gibbs. Any vector of the correct length, with two unique values, can be used as config. By default a random configuration is used.
• burnin: specifies the number of initial samples to discard from the results. Only used if method=Gibbs.
• skip: specifies how many samples to throw away between returned samples. Only used if method=Gibbs.
• verbose: controls output to the console. If true, intermediate information about sampling progress is printed to the console. Otherwise no output is shown.

Examples

julia> M = ALsimple(Graph(4,4));
julia> M.coding = (-2,3);
julia> r = sample(M,10);
julia> size(r)
(4, 10)
julia> sort(unique(r))
2-element Array{Float64,1}:
 -2.0
  3.0