stRing

Here is some more code, to keep the troll confused and happy.


library(ggm)

# Figure 3.1 data
Data=read.table("http://pages.usherbrooke.ca/jshipley/recherche/NINA%20files/fig31.dat")

# Write the DAG
DAG1=DAG(V2~V1, V3~V2, V4~V2, V5~V3, V5~V4, V5~V6)
# Draw the DAG
drawGraph(DAG1)
# D-separation test
dSep(DAG1, "V3","V1", NULL)
# Shipley Test
shipley.test(DAG1, cov(Data), dim(Data)[1])

#####################################################
#### SEM

Data=read.table("http://pages.usherbrooke.ca/jshipley/recherche/NINA%20files/fig31.dat")
names(Data)=LETTERS[1:dim(Data)[2]]

library(sem)
# Input the model
path3.1=specify.model()
A -> B, beta.AB, NA
B -> C, beta.BC, NA
B -> D, beta.BD, NA
C -> E, beta.CE, NA
D -> E, beta.DE, NA
F -> E, beta.FE, NA
A <-> A, var.A, NA
B <-> B, var.B, NA
C <-> C, var.C, NA
D <-> D, var.D, NA
E <-> E, var.E, NA
F <-> F, var.F, NA

# Fit the SEM
model3.1=sem(path3.1, cov(Data), N=dim(Data)[1])
summary(model3.1)

# Specify an alternative model: remove the B->D path
path3.2=specify.model()
A -> B, beta.AB, NA
B -> C, beta.BC, NA
C -> E, beta.CE, NA
D -> E, beta.DE, NA
F -> E, beta.FE, NA
A <-> A, var.A, NA
B <-> B, var.B, NA
C <-> C, var.C, NA
D <-> D, var.D, NA
E <-> E, var.E, NA
F <-> F, var.F, NA

path3.2

model3.2=sem(path3.2, cov(Data), N=dim(Data)[1])
summary(model3.2)

# Compare the two models
anova(model3.1, model3.2)

###################
# String

# Read data (use your own file!)
string=read.table("string.txt")

# First parameterisation
path.string1=specify.model()
Length -> V1, NA, 1
Length -> V2, beta.2, NA
Length -> V3, beta.3, NA
Length -> V4, beta.4, NA
V1 <-> V1, var.1, NA
V2 <-> V2, var.2, NA
V3 <-> V3, var.3, NA
V4 <-> V4, var.4, NA
Length <-> Length, var.L, NA

model.string1=sem(path.string1, cov(string), N=dim(string)[1])
summary(model.string1)

# First parameterisation, Var(V1) set to zero
path.string1a=specify.model()
Length -> V1, NA, 1
Length -> V2, beta.2, NA
Length -> V3, beta.3, NA
Length -> V4, beta.4, NA
V1 <-> V1, NA, 0
V2 <-> V2, var.2, NA
V3 <-> V3, var.3, NA
V4 <-> V4, var.4, NA
Length <-> Length, var.L, NA

model.string1a=sem(path.string1a, cov(string), N=dim(string)[1])
summary(model.string1a)

# Second parameterisation
path.string2=specify.model()
Length -> V1, beta.1, NA
Length -> V2, beta.2, NA
Length -> V3, beta.3, NA
Length -> V4, beta.4, NA
V1 <-> V1, var.1, NA
V2 <-> V2, var.2, NA
V3 <-> V3, var.3, NA
V4 <-> V4, var.4, NA
Length <-> Length, NA, 1

model.string2=sem(path.string2, cov(string), N=dim(string)[1])
summary(model.string2)

# Second parameterisation, Var(V1) set to zero
path.string2a=specify.model()
Length -> V1, beta.1, NA
Length -> V2, beta.2, NA
Length -> V3, beta.3, NA
Length -> V4, beta.4, NA
V1 <-> V1, NA, 0
V2 <-> V2, var.2, NA
V3 <-> V3, var.3, NA
V4 <-> V4, var.4, NA
Length <-> Length, NA, 1

model.string2a=sem(path.string2a, cov(string), N=dim(string)[1])
summary(model.string2a)

# Profile likelihood, to look at likelihood of Var(V1)
GetDev=function(var, path, Cov, N, ...) {
path[5,3] = var
deviance(sem(path, S=Cov, N=N, ...))
}

Var1=as.character(seq(0,5, length=100))
Deviance=sapply(Var1, GetDev, path=path.string1a, Cov=cov(string), N=dim(string)[1])

plot(Var1, -Deviance+min(Deviance), type="l")
abline(h=-3.71, col=2)

DevDiff=-Deviance+min(Deviance)
max(Var1[which(DevDiff>-3.71)])