1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
|
# last modified 2019-11-15 by J. Fox
# Modified for Compiled Objective and nlm in C/C++ by Zhenghua Nie.
optimizerSem <- function(start, objective=objectiveML,
gradient=TRUE, maxiter, debug, par.size, model.description, warn, ...){
with(model.description, {
obj <- objective(gradient=gradient)$objective
typsize <- if (par.size == 'startvalues') abs(start) else rep(1, t) #we move typesize in csem.cpp.
# objectiveCompiled <- "objectiveML"
# if(identical(objective, objectiveCompiledML) || identical(objective, objectiveML))
# objectiveCompiled <- "objectiveML"
# if(identical(objective, objectiveCompiledGLS) || identical(objective, objectiveGLS))
# objectiveCompiled <- "objectiveGLS"
# objectiveCompiled <- deparse(substitute(objective))
# if (!objectiveCompiled %in% c("objectiveML", "objectiveGLS")) stop("optimizerSem requires the objectiveML or objectiveGLS objective function")
# objectiveCompiled <- "objectiveML"
if(identical(objective, objectiveML)) objectiveCompiled <- "objectiveML"
else if (identical(objective, objectiveGLS)) objectiveCompiled <- "objectiveGLS"
else if (identical(objective, objectiveFIML)) objectiveCompiled <- "objectiveFIML"
else stop("optimizerSem requires the objectiveML or objectiveGLS or objectiveFIML objective function")
if (!warn) save.warn <- options(warn=-1)
res <- CompiledSolve(model.description=model.description, start=start, objective=objectiveCompiled, gradient=gradient, typsize=typsize, debug=debug, maxiter=maxiter)
if (!warn) options(save.warn)
result <- list()
result$convergence <- res$code <= 2
result$iterations <- res$iterations
par <- res$estimate
names(par) <- param.names
result$par <- par
if (!result$convergence)
warning(paste('Optimization may not have converged; nlm return code = ',
res$code, '. Consult ?nlm.\n', sep=""))
vcov <- matrix(NA, t, t)
qr.hess <- try(qr(res$hessian), silent=TRUE)
if (inherits(qr.hess, "try-error")){
warning("Could not compute QR decomposition of Hessian.\nOptimization probably did not converge.\n")
}
else if (qr.hess$rank < t){
warning(' singular Hessian: model is probably underidentified.\n')
which.aliased <- qr.hess$pivot[-(1:qr.hess$rank)]
result$aliased <- param.names[which.aliased]
}
else {
vcov <- (2/(N - (!raw))) * solve(res$hessian)
if (any(diag(vcov) < 0)) {
result$aliased <- param.names[diag(vcov) < 0]
warning("Negative parameter variances.\nModel may be underidentified.\n")
}
}
colnames(vcov) <- rownames(vcov) <- param.names
result$vcov <- vcov
result$criterion <- res$minimum # c(result$obj) - n - log(det(S))
# if(identical(objective, objectiveML) || identical(objective, objectiveGLS)){
C <- res$C
A <- res$A
P <- res$P
# } else {
# obj <- obj(par, model.description)
# C <- attr(obj, "C")
# A <- attr(obj, "A")
# P <- attr(obj, "P")
# }
rownames(C) <- colnames(C) <- var.names[observed]
result$C <- C
rownames(A) <- colnames(A) <- var.names
result$A <- A
rownames(P) <- colnames(P) <- var.names
result$P <- P
class(result) <- "semResult"
result
}
)
}
|
