# File R/ergmm.probs.R in package latentnet, part of the
# Statnet suite of packages for network analysis, https://statnet.org .
#
# This software is distributed under the GPL-3 license. It is free,
# open source, and has the attribution requirements (GPL Section 7) at
# https://statnet.org/attribution .
#
# Copyright 2003-2024 Statnet Commons
################################################################################
not.given<-function(name,theta,given){
is.null(given[[name]]) && !is.null(theta[[name]])
}
PRIOR_NAMES<-list(sender=c("sender.var"),
receiver=c("receiver.var"),
sociality=c("sociality.var"),
Z=c("Z.var"))
lp.works<-function(name,theta,given){
not.given(name, theta, given) && all(PRIOR_NAMES[[name]]%in%names(.merge.lists(theta,given)))
}
bipartite.augment<-function(m){
actors<-dim(m)[1]
events<-dim(m)[2]
rbind(matrix(NA,actors,actors+events),
cbind(t(m),matrix(NA,events,events)))
}
getYm<-function(Yg,response=NULL){
m <-
if(is.null(response)){
as.matrix.network(Yg, matrix.type="adjacency")
}else{
if(is.matrix(response)) response
else as.matrix.network(Yg, response, matrix.type="adjacency")
}
## If bipartite, augment into a matrix of the form
## N m
## t(m) N
## where N is a matrix of NAs of appropriate dimension.
if(is.bipartite(Yg)) m<-bipartite.augment(m)
m[!observed.dyads(Yg)]<-NA
m
}
ergmm.eta<-function(model,theta){
n<-network.size(model[["Yg"]])
dir<-is.directed(model[["Yg"]])
eta<-matrix(0,n,n)
if(!is.null(theta[["Z"]]))
eta<-eta+latent.effect.fs[[model[["latentID"]]]](theta[["Z"]])
if(!is.null(theta[["beta"]]))
for(k in seq_along(theta[["beta"]]))
eta<-eta+theta[["beta"]][k]*model[["X"]][[k]]
if(!is.null(theta[["sociality"]])){
eta<-eta+theta[["sociality"]]
eta<-t(t(eta)+theta[["sociality"]])
}
if(!is.null(theta[["sender"]]))
eta<-eta+theta[["sender"]]
if(!is.null(theta[["receiver"]]))
eta<-t(t(eta)+theta[["receiver"]])
return(eta)
}
ergmm.EY<-function(model,theta,NA.unobserved=TRUE){
eta<-ergmm.eta(model,theta)
if(NA.unobserved) eta[!observed.dyads(model[["Yg"]])]<-NA
EY.fs[[model[["familyID"]]]](eta,dispersion=theta[["dispersion"]],fam.par=model[["fam.par"]])
}
ergmm.lpY<-function(model,theta,given=list(),up.to.const=FALSE){
theta<-.merge.lists(theta,given)
Yg<-model[["Yg"]]
Ym<-model[["Ym"]]
n<-network.size(Yg)
eta<-ergmm.eta(model,theta)
obs<-observed.dyads(Yg)
lpY<-if(up.to.const) lpYc.fs[[model[["familyID"]]]](Ym[obs],eta[obs],dispersion=theta[["dispersion"]],model[["fam.par"]]) else lpY.fs[[model[["familyID"]]]](Ym[obs],eta[obs],dispersion=theta[["dispersion"]],model[["fam.par"]])
return(sum(lpY))
}
ergmm.lpY.grad<-function(model,theta,given=list()){
theta<-.merge.lists(theta,given)
n<-network.size(model[["Yg"]])
obs<-observed.dyads(model[["Yg"]])
eta<-ergmm.eta(model,theta)
dlpY.deta <- dlpY.deta.fs[[model[["familyID"]]]](model[["Ym"]],eta,dispersion=theta[["dispersion"]],model[["fam.par"]])
dlpY.deta[!obs] <- 0
grad<-list()
if(not.given("beta",theta,given)) grad[["beta"]] <- sapply(seq_along(theta[["beta"]]),function(k) sum((dlpY.deta*model[["X"]][[k]])[obs]))
if(not.given("Z",theta,given)){
grad[["Z"]]<-dlpY.dZ.fs[[model[["latentID"]]]](theta[["Z"]],dlpY.deta)
}
if(not.given("sociality",theta,given))
grad[["sociality"]] <- sapply(seq_len(n),function(i) sum(dlpY.deta[i,][obs[i,]])+sum(dlpY.deta[,i][obs[,i]]))
else{
if(not.given("sender",theta,given)) grad[["sender"]] <- sapply(seq_len(n),function(i) sum(dlpY.deta[i,][obs[i,]]))
if(not.given("receiver",theta,given)) grad[["receiver"]] <- sapply(seq_len(n),function(i) sum(dlpY.deta[,i][obs[,i]]))
}
if(not.given("dispersion",theta,given)){
grad[["dispersion"]] <- dlpY.ddispersion.fs[[model[["familyID"]]]](model[["Ym"]],eta,dispersion=theta[["dispersion"]],model[["fam.par"]])
}
grad
}
ergmm.lpY.C<-function(model,theta){
Y <- model[["Ym"]]
Y[is.na(Y)] <- 0
n <- network.size(model[["Yg"]])
## Figure out the design matrix.
observed<-observed.dyads(model[["Yg"]])
if((observed==(diag(n)==0) && is.directed(model[["Yg"]])) ||
(observed==lower.tri(diag(n)) && !is.directed(model[["Yg"]])))
observed<-NULL
## Sanity checks: the following block of code checks that all dimensionalities and
## dimensions are correct, and those optional parameters that are required by the presence
## of other optional parameters are present.
for(i in seq_len(model[["p"]]))
if(!all(dim(model[["X"]][[i]])==c(n,n))) stop("Incorrect size for covariate matrices.")
if(!is.null(theta[["Z"]])){
if(!all(dim(theta[["Z"]])==c(n,model[["d"]]))) stop("Incorrect size for the latent positions.")
}
if(length(theta[["beta"]])!=model[["p"]]) stop("Incorrect length for the beta vector.")
if(!is.null(theta[["sociality"]])){
if(length(theta[["sociality"]])!=n) stop("Incorrect length for the vector of sociality effects.")
}
if(!is.null(theta[["sender"]])){
if(length(theta[["sender"]])!=n) stop("Incorrect length for the vector of sender effects.")
}
if(!is.null(theta[["receiver"]])){
if(length(theta[["receiver"]])!=n) stop("Incorrect length for the vector of receiver effects.")
}
## End Sanity checks.
ret <- .C("ERGMM_lp_Y_wrapper",
n=as.integer(n), p=as.integer(model[["p"]]),
d=as.integer(model[["d"]]), latent=as.integer(NVL(model[["latentID"]],0)), family=as.integer(NVL(model[["familyID"]],0)), res=as.integer(with(model,c(sender,receiver,sociality))),
dir=as.integer(is.directed(model[["Yg"]])),
viY=as.integer(Y),
vdY=as.double(Y),
iconsts=as.integer(model[["iconsts"]]), dconsts=as.double(model[["dconsts"]]),
vX=as.double(unlist(model[["X"]])),
Z=as.double(theta[["Z"]]),
beta=as.double(theta[["beta"]]),
sender=if(is.null(theta[["sociality"]])) as.double(theta[["sender"]]) else as.double(theta[["sociality"]]), receiver=as.double(theta[["receiver"]]),
dispersion=NVL(theta[["dispersion"]],0),
observed=as.integer(NVL(observed,-1)),
lpY=double(1),
PACKAGE="latentnet")
ret[["lpY"]]
}
observed.dyads<-function(Yg){
observed.dyads<-get.network.attribute(Yg,"design")
if(is.null(observed.dyads)){
if(!is.bipartite(Yg))
observed.dyads<-!is.na(as.matrix.network(Yg,matrix.type="adjacency"))
else
observed.dyads<-bipartite.augment(!is.na(as.matrix.network(Yg,matrix.type="adjacency")))
}
else{
observed.dyads<-bipartite.augment(as.matrix.network(observed.dyads,matrix.type="adjacency")==0)
}
observed.dyads[is.na(observed.dyads)]<-FALSE
if(!is.directed(Yg)) observed.dyads[upper.tri(observed.dyads)]<-FALSE
if(!is.bipartite(Yg) && !has.loops(Yg)) diag(observed.dyads)<-FALSE
observed.dyads
}
pack.optim<-function(theta,fit.vars=NULL){
if(is.null(fit.vars))
return(c(theta[["beta"]],theta[["Z"]],
theta[["sender"]],theta[["receiver"]],theta[["sociality"]],
theta[["Z.var"]],theta[["Z.mean"]],
theta[["sender.var"]],theta[["receiver.var"]],theta[["sociality.var"]],
theta[["dispersion"]]))
else
return(c(if(fit.vars[["beta"]])theta[["beta"]],
if(fit.vars[["Z"]])theta[["Z"]],
if(fit.vars[["sender"]])theta[["sender"]],
if(fit.vars[["receiver"]])theta[["receiver"]],
if(fit.vars[["sociality"]])theta[["sociality"]],
if(fit.vars[["Z.var"]])theta[["Z.var"]],
if(fit.vars[["Z.mean"]])theta[["Z.mean"]],
if(fit.vars[["sender.var"]])theta[["sender.var"]],
if(fit.vars[["receiver.var"]])theta[["receiver.var"]],
if(fit.vars[["sociality.var"]])theta[["sociality.var"]],
if(fit.vars[["dispersion"]])theta[["dispersion"]]))
}
reg.fit.vars<-function(fit.vars){
for(name in ERGMM.PAR_VAR_NAMES)
if(!(name %in% names(fit.vars)))fit.vars[[name]]<-FALSE
fit.vars
}
inv.fit.vars<-function(fit.vars){
for(name in ERGMM.PAR_VAR_NAMES)
if(name %in% names(fit.vars))
fit.vars[[name]]<-!fit.vars[[name]]
fit.vars
}
FIT_ALL<-list(beta=TRUE,Z=TRUE,sender=TRUE,receiver=TRUE,sociality=TRUE,
Z.var=TRUE,Z.mean=TRUE,
sender.var=TRUE,receiver.var=TRUE,sociality.var=TRUE,dispersion=TRUE)
FIT_MLE<-list(beta=TRUE,Z=TRUE,sender=TRUE,receiver=TRUE,sociality=TRUE,dispersion=TRUE)
unpack.optim<-function(v,fit.vars,model){
p<-model[["p"]]
n<-network.size(model[["Yg"]])
G<-model[["G"]]
d<-model[["d"]]
v.must.be<-with(fit.vars,
beta*p +
Z*n*d +
sender*n +
receiver*n +
sociality*n +
Z.var*(d>0)*max(1,G) +
Z.mean*G*d +
sender.var +
receiver.var +
sociality.var +
dispersion)
if(length(v)!=v.must.be){
stop(paste("Input vector wrong length: ", length(v),
" but should be ",v.must.be,".",
sep=""))
}
pos<-0
ret<-list()
if(fit.vars[["beta"]] && p>0){
ret[["beta"]]<-v[pos+seq_len(p)]
pos<-pos+p
}
if(fit.vars[["Z"]] && d>0){
ret[["Z"]]<-matrix(v[pos+seq_len(n*d)],nrow=n,ncol=d)
pos<-pos+n*d
}
if(fit.vars[["sender"]]){
ret[["sender"]]<-v[pos+seq_len(n)]
pos<-pos+n
}
if(fit.vars[["receiver"]]){
ret[["receiver"]]<-v[pos+seq_len(n)]
pos<-pos+n
}
if(fit.vars[["sociality"]]){
ret[["sociality"]]<-v[pos+seq_len(n)]
pos<-pos+n
}
if(fit.vars[["Z.var"]] && d>0){
ret[["Z.var"]]<-v[pos+seq_len(max(1,G))]
pos<-pos+max(1,G)
}
if(fit.vars[["Z.mean"]] && d>0 && G>0){
ret[["Z.mean"]]<-matrix(v[pos+seq_len(G*d)],nrow=G,ncol=d)
pos<-pos+G*d
}
if(fit.vars[["sender.var"]]){
ret[["sender.var"]]<-v[pos+1]
pos<-pos+1
}
if(fit.vars[["receiver.var"]]){
ret[["receiver.var"]]<-v[pos+1]
pos<-pos+1
}
if(fit.vars[["sociality.var"]]){
ret[["sociality.var"]]<-v[pos+1]
pos<-pos+1
}
if(fit.vars[["dispersion"]]){
ret[["dispersion"]]<-v[pos+1]
pos<-pos+1
}
ret
}
mk.lp.optim.fs<-function(fit.vars,model,prior,given=list(),opt=c("lpY","lpZ","lpBeta","lpRE","lpREV","lpLV","lpdispersion")){
fit.vars<-reg.fit.vars(fit.vars)
return(list(
f=function(v){
theta<-unpack.optim(v,fit.vars,model)
ergmm.lp(model,theta,prior=prior,given=given,
opt=opt,
up.to.const=TRUE)
},
grad.f=function(v){
theta<-unpack.optim(v,fit.vars,model)
# note that ergmm.lp.grad doesn't take the up.to.const parameter
gr<-ergmm.lp.grad(model,theta,prior=prior,given=given,
opt=opt)
pack.optim(gr,fit.vars)
}
)
)
}
find.mle<-function(model,start,given=list(),control,
hessian=FALSE,mllk=TRUE){
fit.vars<-list()
for(name in ERGMM.PAR_LLK_NAMES)
fit.vars[[name]]<-not.given(name,start,given)
mpe<-find.mpe(model,start,given=given,control=control,
hessian=hessian,mlp=mllk,opt="lpY",fit.vars=fit.vars)
if(mllk) mpe[["lpY"]]<-mpe[["mlp"]]
mpe
}
find.mpe<-function(model,start,given=list(),prior=list(),control,fit.vars=NULL,opt=c("lpY","lpZ","lpBeta","lpRE","lpREV","lpLV","lpdispersion"),
hessian=FALSE,mlp=TRUE){
if(is.null(fit.vars)){
fit.vars<-list()
for(name in names(start))
fit.vars[[name]]<-not.given(name,start,given)
}else{
fit.vars<-reg.fit.vars(fit.vars)
for(name in names(fit.vars))
if(!not.given(name,start,given)) fit.vars[[name]]<-FALSE
}
fit.vars<-reg.fit.vars(fit.vars)
optim.control<-list(fnscale=-1,
maxit=control[["mle.maxit"]],
trace=max(0,control[["verbose"]]-2))
optim.fs<-mk.lp.optim.fs(fit.vars,model,prior=prior,given=given,opt=opt)
start.vals<-pack.optim(start,fit.vars)
p<-model[["p"]]
n<-network.size(model[["Yg"]])
G<-model[["G"]]
d<-model[["d"]]
#' @importFrom stats optim
vmpe <- ##try(
optim(par=start.vals,fn=optim.fs[["f"]],gr=optim.fs[["grad.f"]],
method="L-BFGS-B",
lower=pack.optim(list(
beta=rep(-Inf,p),
Z=rep(-Inf,n*d),
sender=rep(-Inf,n),
receiver=rep(-Inf,n),
sociality=rep(-Inf,n),
Z.var=rep(sqrt(.Machine[["double.eps"]]),(d>0)*max(G,1)),
Z.mean=rep(-Inf,d*G),
sender.var=sqrt(.Machine[["double.eps"]]),
receiver.var=sqrt(.Machine[["double.eps"]]),
sociality.var=sqrt(.Machine[["double.eps"]]),
dispersion=sqrt(.Machine[["double.eps"]])),
fit.vars=fit.vars),
control=optim.control,hessian=hessian)
## )
if(inherits(vmpe,"try-error")) return(NULL)
mpe<-unpack.optim(vmpe[["par"]],fit.vars,model)
mpe<-.merge.lists(mpe,given)
mpe[["Z.K"]]<-.merge.lists(start,given)[["Z.K"]]
mpe[["Z.pK"]]<-if(!is.null(mpe[["Z.K"]])) tabulate(mpe[["Z.K"]],G)/n
if(mlp)
mpe[["mlp"]]<-ergmm.lp(model,mpe,prior=prior,given=given,opt=opt)
if(hessian) mpe[["hessian"]]<-vmpe[["hessian"]]
mpe
}
ergmm.lp<-function(model,theta,prior,given=list(),opt=c("lpY","lpZ","lpBeta","lpRE","lpREV","lpLV"),up.to.const=FALSE){
lpY<-if("lpY" %in% opt) ergmm.lpY(model,theta,
given=given,up.to.const=up.to.const) else 0
lpZ<-if("lpZ" %in% opt) ergmm.lpZ(theta,given=given) else 0
lpRE<-if("lpRE" %in% opt) ergmm.lpRE(theta,given=given) else 0
lpBeta<-if("lpBeta" %in% opt) ergmm.lpBeta(theta,prior,given=given) else 0
lpREV<-if("lpREV" %in% opt) ergmm.lpREV(theta,prior,given=given) else 0
lpLV<-if("lpLV" %in% opt) ergmm.lpLV(theta,prior,given=given) else 0
lpdispersion<-if("lpdispersion" %in% opt) ergmm.lpdispersion(theta,prior,given=given) else 0
lpAll<-lpY+lpZ+lpRE+lpBeta+lpREV+lpLV+lpdispersion
return(lpAll)
}
.merge.lists<-function(...){
out<-list()
for(l in list(...)){
for(name in names(l))
out[[name]]<-l[[name]]
if(!inherits(l,"list"))
class(out)<-class(l)
}
out
}
.sum.lists<-function(...){
out<-list()
for(l in list(...)){
for(name in names(l)){
if(name %in% names(out))
out[[name]]<-out[[name]]+l[[name]]
else
out[[name]]<-l[[name]]
}
}
class(out)<-class(list(l))
out
}
zero.list<-function(x){
out<-x
for(name in names(x)){
d<-dim(x[[name]])
if(is.null(d)) out[[name]]<-rep(0,length(x[[name]]))
else out[[name]]<-array(0,dim=d)
}
}
filter.list<-function(x,y){
out<-x
for(name in names(y)){
if(!isTRUE(y)) out[[name]]<-NULL
}
out
}
cmp.lists<-function(x,y){
out<-list()
for(name in names(x))
if(name %in% names(y))
out[[name]]<-mean(abs(x[[name]]-y[[name]]))/mean(abs(x[[name]]))
out
}
ergmm.lp.grad<-function(model,theta,prior,given=list(),opt=c("lpY","lpZ","lpBeta","lpRE","lpREV","lpLV","lpdispersion")){
grad<-.sum.lists(if("lpY" %in% opt) if(not.given("beta",theta,given)||
not.given("Z",theta,given)||
not.given("sender",theta,given)||
not.given("receiver",theta,given)||
not.given("sociality",theta,given)||
not.given("dispersion",theta,given)) ergmm.lpY.grad(model,theta,given=given),
if("lpZ" %in% opt) ergmm.lpZ.grad(theta,given=given),
if("lpRE" %in% opt) ergmm.lpRE.grad(theta,given=given),
if("lpBeta" %in% opt) ergmm.lpBeta.grad(theta,prior,given=given),
if("lpREV" %in% opt) ergmm.lpREV.grad(theta,prior,given=given),
if("lpLV" %in% opt) ergmm.lpLV.grad(theta,prior,given=given),
if("lpdispersion" %in% opt) ergmm.lpdispersion.grad(theta,prior,given=given))
grad
}
ergmm.lp.grad.approx<-function(which.vars,model,theta,prior,delta,given=list(),opt=c("lpY","lpZ","lpBeta","lpRE","lpREV","lpLV","lpdispersion")){
which.vars[["Z.K"]]<-FALSE
which.vars<-reg.fit.vars(which.vars)
for(var in names(which.vars)) if(!(var %in% names(theta))) which.vars[[var]]<-FALSE
v<-pack.optim(theta,which.vars)
dlpdv<-numeric(length(v))
for(i in seq_along(v)){
v.m<-v.p<-v
v.m[i]<-v.m[i]-delta
theta.m<-.merge.lists(theta,unpack.optim(v.m,which.vars,model))
lp.m<-ergmm.lp(model,theta.m,prior,given=given,opt=opt)
v.p[i]<-v.p[i]+delta
theta.p<-.merge.lists(theta,unpack.optim(v.p,which.vars,model))
lp.p<-ergmm.lp(model,theta.p,prior,given=given,opt=opt)
dlpdv[i]<-(lp.p-lp.m)/(2*delta)
}
return(unpack.optim(dlpdv,which.vars,model))
}
ergmm.lpZ<-function(theta,given=list()){
theta<-.merge.lists(theta,given)
if(lp.works("Z",theta,given)){
n<-dim(theta[["Z"]])[1]
d<-dim(theta[["Z"]])[2]
if(is.null(theta[["Z.K"]])){
if(!is.null(theta[["Z.mean"]])) stop("Given cluster means without cluster assignments!")
theta[["Z.K"]]<-rep(1,n)
theta[["Z.mean"]]<-matrix(0,nrow=1,ncol=d)
}
sum(dnorm(theta[["Z"]],theta[["Z.mean"]][theta[["Z.K"]],],matrix(sqrt(theta[["Z.var"]][theta[["Z.K"]]]),nrow=n,ncol=d,byrow=FALSE),TRUE))
}
else 0
}
ergmm.lpZ.grad<-function(theta,given=list()){
theta<-.merge.lists(theta,given)
deriv<-list()
if(lp.works("Z",theta,given)){
n<-dim(theta[["Z"]])[1]
d<-dim(theta[["Z"]])[2]
G<-if(is.null(theta[["Z.K"]])) 0 else dim(theta[["Z.mean"]])[1]
if(is.null(theta[["Z.K"]])){
if(!is.null(theta[["Z.mean"]])) stop("Given cluster means without cluster assignments!")
theta[["Z.K"]]<-rep(1,n)
theta[["Z.mean"]]<-matrix(0,nrow=1,ncol=d)
}
Z.dev<-(theta[["Z"]]-theta[["Z.mean"]][theta[["Z.K"]],])/matrix(theta[["Z.var"]][theta[["Z.K"]]],nrow=n,ncol=d,byrow=FALSE)
deriv[["Z"]]<--Z.dev
if(not.given("Z.var",theta,given)) deriv[["Z.var"]]<-sapply(seq_len(max(G,1)),
function(g)
(sum(Z.dev[theta[["Z.K"]]==g,,drop=FALSE]^2)-d*sum(theta[["Z.K"]]==g)/theta[["Z.var"]][g])/2)
if(not.given("Z.mean",theta,given) && G) deriv[["Z.mean"]]<-do.call(rbind,lapply(seq_len(G),function(g) apply(Z.dev[theta[["Z.K"]]==g,,drop=FALSE],2,sum)))
}
deriv
}
ergmm.lpRE<-function(theta, given=list()){
theta<-.merge.lists(theta,given)
({if(lp.works("sender",theta,given)) sum(dnorm(theta[["sender"]],0,sqrt(theta[["sender.var"]]),TRUE)) else 0}+
{if(lp.works("receiver",theta,given)) sum(dnorm(theta[["receiver"]],0,sqrt(theta[["receiver.var"]]),TRUE)) else 0}+
{if(lp.works("sociality",theta,given)) sum(dnorm(theta[["sociality"]],0,sqrt(theta[["sociality.var"]]),TRUE)) else 0})
}
ergmm.lpRE.grad<-function(theta,given=list()){
theta<-.merge.lists(theta,given)
deriv<-list()
if(lp.works("sender",theta,given)) deriv[["sender"]]<--theta[["sender"]]/theta[["sender.var"]]
if(not.given("sender.var",theta,given)) deriv[["sender.var"]]<-(sum(theta[["sender"]]^2)/theta[["sender.var"]]-length(theta[["sender"]]))/theta[["sender.var"]]/2
if(lp.works("receiver",theta,given)) deriv[["receiver"]]<--theta[["receiver"]]/theta[["receiver.var"]]
if(not.given("receiver.var",theta,given)) deriv[["receiver.var"]]<-(sum(theta[["receiver"]]^2)/theta[["receiver.var"]]-length(theta[["receiver"]]))/theta[["receiver.var"]]/2
if(lp.works("sociality",theta,given)) deriv[["sociality"]]<--theta[["sociality"]]/theta[["sociality.var"]]
if(not.given("sociality.var",theta,given)) deriv[["sociality.var"]]<-(sum(theta[["sociality"]]^2)/theta[["sociality.var"]]-length(theta[["sociality"]]))/theta[["sociality.var"]]/2
deriv
}
ergmm.lpBeta<-function(theta,prior,given=list()){
theta<-.merge.lists(theta,given)
return({if(not.given("beta",theta,given)) sum(dnorm(theta[["beta"]],prior[["beta.mean"]],sqrt(prior[["beta.var"]]),TRUE)) else 0})
}
ergmm.lpBeta.grad<-function(theta,prior,given=list()){
theta<-.merge.lists(theta,given)
deriv<-list()
if(not.given("beta",theta,given)) deriv[["beta"]]<--(theta[["beta"]]-prior[["beta.mean"]])/prior[["beta.var"]]
deriv
}
#' @importFrom stats dchisq
dsclinvchisq<-function(x,df,scale=1,log=FALSE){
if(log) dchisq(df*scale/x,df,log=TRUE)+log(df)+log(scale)-2*log(x)
else dchisq(df*scale/x,df,log=FALSE)*df*scale/x/x
}
ergmm.lpREV<-function(theta,prior,given=list()){
theta<-.merge.lists(theta,given)
({if(not.given("sender.var",theta,given)) dsclinvchisq(theta[["sender.var"]],prior[["sender.var.df"]],prior[["sender.var"]],TRUE) else 0}+
{if(not.given("receiver.var",theta,given)) dsclinvchisq(theta[["receiver.var"]],prior[["receiver.var.df"]],prior[["receiver.var"]],TRUE) else 0}+
{if(not.given("sociality.var",theta,given)) dsclinvchisq(theta[["sociality.var"]],prior[["sociality.var.df"]],prior[["sociality.var"]],TRUE) else 0})
}
ergmm.lpREV.grad<-function(theta,prior,given=list()){
theta<-.merge.lists(theta,given)
deriv<-list()
if(not.given("sender.var",theta,given)) deriv[["sender.var"]]<-prior[["sender.var.df"]]*prior[["sender.var"]]/theta[["sender.var"]]^2/2-(prior[["sender.var.df"]]/2+1)/theta[["sender.var"]]
if(not.given("receiver.var",theta,given)) deriv[["receiver.var"]]<-prior[["receiver.var.df"]]*prior[["receiver.var"]]/theta[["receiver.var"]]^2/2-(prior[["receiver.var.df"]]/2+1)/theta[["receiver.var"]]
if(not.given("sociality.var",theta,given)) deriv[["sociality.var"]]<-prior[["sociality.var.df"]]*prior[["sociality.var"]]/theta[["sociality.var"]]^2/2-(prior[["sociality.var.df"]]/2+1)/theta[["sociality.var"]]
deriv
}
ergmm.lpLV<-function(theta,prior,given=list()){
theta<-.merge.lists(theta,given)
({if(not.given("Z.var",theta,given)) sum(dsclinvchisq(theta[["Z.var"]],prior[["Z.var.df"]],prior[["Z.var"]],log=TRUE)) else 0}+
{if(not.given("Z.mean",theta,given)) sum(dnorm(theta[["Z.mean"]],0,sqrt(prior[["Z.mean.var"]]),log=TRUE)) else 0})
}
ergmm.lpLV.grad<-function(theta,prior,given=list()){
theta<-.merge.lists(theta,given)
deriv<-list()
if(not.given("Z.var",theta,given)) deriv[["Z.var"]]<-prior[["Z.var.df"]]*prior[["Z.var"]]/theta[["Z.var"]]^2/2-(prior[["Z.var.df"]]/2+1)/theta[["Z.var"]]
if(not.given("Z.mean",theta,given)) deriv[["Z.mean"]]<--theta[["Z.mean"]]/prior[["Z.mean.var"]]
deriv
}
ergmm.lpdispersion<-function(theta,prior,given=list()){
theta<-.merge.lists(theta,given)
(if(not.given("dispersion",theta,given)) dsclinvchisq(theta[["dispersion"]],prior[["dispersion.df"]],prior[["dispersion"]],TRUE) else 0)
}
ergmm.lpdispersion.grad<-function(theta,prior,given=list()){
theta<-.merge.lists(theta,given)
deriv<-list()
if(not.given("dispersion",theta,given)) deriv[["dispersion"]]<-prior[["dispersion.df"]]*prior[["dispersion"]]/theta[["dispersion"]]^2/2-(prior[["dispersion.df"]]/2+1)/theta[["dispersion"]]
deriv
}
lpsum<-function(theta,which=c("lpY","lpZ","lpBeta","lpRE","lpREV","lpLV","lpdispersion")){
sum(sapply(which,function(lp) if(is.null(theta[[lp]])) 0 else theta[[lp]]))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.