tests/regression_tests/hojsgaard_model_tests/triangle_model/triangle-model_loglin-gRim.R

In npetraco/CRFutil: Handy Utility Functions for Use with CRF and gRbase packages

# loglin gRim
library(gRim)
library(MASS)

# Load samps generated by true model:
fpth <- "/home/npetraco/codes/R/CRFutil/tests/regression_tests/hojsgaard_model_tests/triangle_model/triangle_data/"
load(paste0(fpth,"triangle_samp.RData"))
head(samps)

# Fold samples into contingency table for a look
X.cont <- xtabs(~., data=data.frame(samps))
X.cont

# Triangle model
colnames(samps)
grphf <- ~1:2+1:3+2:3
# Check the graph:
gp <- ug(grphf, result = "graph")
dev.off()
iplot(gp)


# Fit the Hojgaard model with loglm and gRim
fact.grphf <- ~X.1 + X.2 + X.3 + X.1:X.2 + X.1:X.3 + X.2:X.3
loglin.mle <- loglm(fact.grphf, data=X.cont); # loglm from MASS which uses loglin in base
loglin.mle

X.loglm.coefs <- coef(loglin.mle)
X.loglm.coefs

# Fit contingency table:
X.cont.fitted <- fitted(loglin.mle)

# Flatten fit contingency table into matrix-table form
X.counts.fitted <- data.frame(ftable(X.cont.fitted))
X.counts.fitted
sum(X.counts.fitted[,4]) # Equal to the sample size?

# Estimate state probabilities is by the loglin fitted relative frquencies:
X.freq.fitted <- X.counts.fitted
X.freq.fitted[,4] <- X.freq.fitted[,4]/sum(X.counts.fitted[,4])
X.freq.fitted <- cbind(X.freq.fitted, X.counts.fitted[,4]) # Carry along the fit counts the computation is based on
colnames(X.freq.fitted)[c(4,5)] <- c("LogLin.Freq","LogLin.Counts")
loglin.dist.info <- X.freq.fitted
save(loglin.dist.info, file = paste0(fpth,"triangle_loglin_dist.RData"))