Last updated on 2020-02-19 10:48:56 CET.
| Flavor | Version | Tinstall | Tcheck | Ttotal | Status | Flags |
|---|---|---|---|---|---|---|
| r-devel-linux-x86_64-debian-clang | 0.4 | 36.76 | 73.74 | 110.50 | ERROR | |
| r-devel-linux-x86_64-debian-gcc | 0.4 | 27.30 | 57.56 | 84.86 | ERROR | |
| r-devel-linux-x86_64-fedora-clang | 0.4 | 141.89 | ERROR | |||
| r-devel-linux-x86_64-fedora-gcc | 0.4 | 137.34 | ERROR | |||
| r-devel-windows-ix86+x86_64 | 0.4 | 68.00 | 121.00 | 189.00 | NOTE | |
| r-devel-windows-ix86+x86_64-gcc8 | 0.4 | 101.00 | 156.00 | 257.00 | ERROR | |
| r-patched-linux-x86_64 | 0.4 | 30.91 | 67.21 | 98.12 | NOTE | |
| r-patched-solaris-x86 | 0.4 | 170.60 | NOTE | |||
| r-release-linux-x86_64 | 0.4 | 28.60 | 66.59 | 95.19 | NOTE | |
| r-release-windows-ix86+x86_64 | 0.4 | 73.00 | 116.00 | 189.00 | NOTE | |
| r-release-osx-x86_64 | 0.4 | NOTE | ||||
| r-oldrel-windows-ix86+x86_64 | 0.4 | 54.00 | 120.00 | 174.00 | NOTE | |
| r-oldrel-osx-x86_64 | 0.4 | NOTE |
Version: 0.4
Check: DESCRIPTION meta-information
Result: NOTE
Malformed Title field: should not end in a period.
Flavors: r-devel-linux-x86_64-debian-clang, r-devel-linux-x86_64-debian-gcc, r-devel-linux-x86_64-fedora-clang, r-devel-linux-x86_64-fedora-gcc, r-devel-windows-ix86+x86_64, r-devel-windows-ix86+x86_64-gcc8, r-patched-linux-x86_64, r-patched-solaris-x86, r-release-linux-x86_64, r-release-windows-ix86+x86_64, r-release-osx-x86_64, r-oldrel-windows-ix86+x86_64, r-oldrel-osx-x86_64
Version: 0.4
Check: R code for possible problems
Result: NOTE
GenerateMixData: no visible global function definition for 'rnbinom'
GenerateMixData: no visible global function definition for 'rbeta'
GenerateMixData: no visible global function definition for 'rbinom'
GenerateMixData: no visible global function definition for 'rpois'
InitializePars: no visible global function definition for 'aggregate'
InitializePars: no visible global function definition for 'rnorm'
MultiStart: no visible global function definition for 'rnorm'
MultiStart: no visible global function definition for 'runif'
inarmix: no visible global function definition for 'terms'
inarmix: no visible global function definition for 'model.frame'
inarmix: no visible global function definition for 'model.response'
inarmix: no visible global function definition for 'model.matrix'
print.diagnose.inarmix: no visible global function definition for
'plot'
print.diagnose.inarmix: no visible global function definition for
'lines'
print.summary.inarmix: no visible global function definition for
'printCoefmat'
summary.inarmix: no visible global function definition for 'pnorm'
Undefined global functions or variables:
aggregate lines model.frame model.matrix model.response plot pnorm
printCoefmat rbeta rbinom rnbinom rnorm rpois runif terms
Consider adding
importFrom("graphics", "lines", "plot")
importFrom("stats", "aggregate", "model.frame", "model.matrix",
"model.response", "pnorm", "printCoefmat", "rbeta",
"rbinom", "rnbinom", "rnorm", "rpois", "runif", "terms")
to your NAMESPACE file.
Flavors: r-devel-linux-x86_64-debian-clang, r-devel-linux-x86_64-debian-gcc, r-devel-linux-x86_64-fedora-clang, r-devel-linux-x86_64-fedora-gcc, r-devel-windows-ix86+x86_64, r-devel-windows-ix86+x86_64-gcc8, r-patched-linux-x86_64, r-patched-solaris-x86, r-release-linux-x86_64, r-release-windows-ix86+x86_64, r-release-osx-x86_64, r-oldrel-windows-ix86+x86_64, r-oldrel-osx-x86_64
Version: 0.4
Check: examples
Result: ERROR
Running examples in 'inarmix-Ex.R' failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: inarmix
> ### Title: Finite mixture model for longitudinal count data.
> ### Aliases: inarmix
> ### Keywords: models regression
>
> ### ** Examples
>
> set.seed(4297)
>
> ############################################################
> #### Simulate data from a two class model
>
> XX <- cbind(rep(1,9),c(0:8)/4)
> colnames(XX) <- c("const","time")
> beta <- rbind(c(-.2,0),c(1.2,.3))
> ### this means that for group 1: (beta_{0},beta_{1}) = (-.2,0)
> ### and for group 2: (beta_{0},beta_{1}) = (1.2,.3)
> autocorr <- c(.2,.2)
> scale <- c(2,2)
> mix.prop <- c(.8,.2) ## proportion in group 1 is .8
>
> testdat <- GenerateMixData(500,beta,autocorr,scale,mix.prop,XX)
> testdat[1:5,]
y subject const time
1 1 1 1 0.00
2 0 1 1 0.25
3 2 1 1 0.50
4 4 1 1 0.75
5 0 1 1 1.00
>
>
> ########################################################################
> #### Fit a linear curve with two classes (with a maximum of 4 iterations)
>
> twoclassfit <- inarmix(y~time,nclasses=2,id=subject,data=testdat,maxiter=4)
Iteration: 1
Iteration: 2
Iteration: 3
Iteration: 4
> summary(twoclassfit)
Call:
inarmix(formula = y ~ time, nclasses = 2, id = subject, data = testdat,
maxiter = 4)
Class Probabilities
Estimate Std.Err
Group1 0.7977567 0.01844
Group2 0.2022433 NA
------------------------------
Group 1 coefficients:
Estimate Std.Err Z value Pr(>z)
(Intercept) -0.18180172 0.05623476 -3.23291 0.0012254 **
time -0.02653325 0.04532190 -0.58544 0.5582519
autocorr. 0.20233773 0.02368218 8.54388 < 2.22e-16 ***
scale 2.00629134 0.08277197 12.15739 < 2.22e-16 ***
------------------------------
Group 2 coefficients:
Estimate Std.Err Z value Pr(>z)
(Intercept) 1.21738430 0.05829181 20.88431 < 2.22e-16 ***
time 0.28683640 0.04407928 6.50728 7.6522e-11 ***
autocorr. 0.19696018 0.03623996 5.43489 5.4830e-08 ***
scale 2.12021748 0.11447226 9.78593 < 2.22e-16 ***
------------------------------
log-likelihood: -6823.222
BIC: 13722.151
AIC: 13664.444
Number of EM iterations: 4
>
> diagnose(twoclassfit)
Did not converge!!!!
Convergence of GEE equations for
each iteration: 1 - converged, 0 - did not converge
iteration GEE 1 GEE 2 l[i] l[i+1]-l[i] ||Psi||^2
[1,] 0 1 1 -6866.134 0.000000 0.646708992
[2,] 1 1 1 -6833.577 -32.556835 0.108210548
[3,] 2 1 1 -6824.850 -8.726774 0.015731280
[4,] 3 1 1 -6823.222 -1.627973 0.002018329
>
> #############################################################
> ##### Fit the same model with specified starting values.
>
> inpars <- list()
> inpars$coef <- rbind(c(-.5,.1),c(.5,0))
> inpars$autocorr <- rep(.3,2)
> inpars$scale <- rep(2,2)
> inpars$mix.prop <- c(.6,.4)
>
> twoclassfit2 <- inarmix(y~time,nclasses=2,id=subject,data=testdat,initvals=inpars,
+ maxiter=4)
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
inarmix
--- call from context ---
inarmix(y ~ time, nclasses = 2, id = subject, data = testdat,
initvals = inpars, maxiter = 4)
--- call from argument ---
if (!check.inits) {
stop("Initial values are not in the correct format")
}
--- R stacktrace ---
where 1: inarmix(y ~ time, nclasses = 2, id = subject, data = testdat,
initvals = inpars, maxiter = 4)
--- value of length: 2 type: logical ---
[1] FALSE TRUE
--- function from context ---
function (formula, nclasses = 1, id, data, initvals = NULL, maxiter = 200,
stoptol = 1e-05, num.restarts = 1, time = NULL, dropthresh = 0.01)
{
mod <- match.call(expand.dots = FALSE)
modterms <- terms(formula)
subj.col <- which(colnames(data) == mod$id)
if (!is.null(mod$time)) {
time.col <- which(colnames(data) == mod$time)
if (length(time.col) == 0) {
stop("time variable not found")
}
}
if (length(subj.col) == 0) {
stop("id variable not found")
}
if (any(is.na(data))) {
stop("There are missing data points")
}
if (is.null(mod$time)) {
data <- data[order(data[, subj.col]), ]
}
else {
data <- data[order(data[, subj.col], data[, time.col]),
]
}
modframe <- model.frame(formula, data)
yy <- model.response(modframe)
XX <- model.matrix(formula, modframe)
QR <- qr(XX)
if (QR$rank < ncol(XX)) {
stop("The design matrix is rank-deficient")
}
if (any(is.infinite(yy)) | any(is.infinite(XX))) {
}
subject.id <- data[, subj.col]
m <- length(unique(subject.id))
nvec <- as.vector(table(subject.id))
loglik.record <- numeric(num.restarts)
initvals.record <- final.record <- list()
ylist <- Xlist <- list()
subject.names <- rep(" ", m)
ncum <- 1 + c(0, cumsum(nvec))
for (i in 1:m) {
ind1 <- ncum[i]
ind2 <- ncum[i + 1] - 1
ylist[[i]] <- yy[ind1:ind2]
Xlist[[i]] <- as.matrix(XX[ind1:ind2, ])
tmp <- unique(data[ind1, subj.col])
subject.names[i] <- as.character(tmp)
}
if (nclasses == 1) {
oneclass.ests <- OneclassEsts(yy, XX, nvec, m)
covs <- OuterProd(beta = oneclass.ests$beta, alpha = oneclass.ests$alpha,
gamma = oneclass.ests$gamma, nvec = nvec, ylist = ylist,
Xlist = Xlist)
cov.mat <- tcrossprod(solve(covs$D, covs$V), solve(covs$D))
std.errs <- sqrt(diag(cov.mat))
names(std.errs) <- c(colnames(XX), "autocorr.", "scale")
betalist <- list()
betalist[[1]] <- oneclass.ests$beta
tmp.pp <- PostProbs(betalist, oneclass.ests$alpha, oneclass.ests$gamma,
1, ylist, Xlist, m, 1)
post.probs <- tmp.pp[[1]]
loglikhood <- tmp.pp[[2]]
results <- list()
class(results) <- "inarmix"
results$call <- mod
results$coefficients <- oneclass.ests$beta
results$alpha <- oneclass.ests$alpha
results$gamma <- oneclass.ests$gamma
results$mix.prop <- NULL
results$coefnames <- colnames(XX)
results$nclasses <- 1
results$post.probs <- NULL
results$niter <- NULL
results$GEE.conv <- NULL
results$loglikhood <- loglikhood
results$em.converged <- NULL
results$pss <- NULL
results$initvals <- NULL
results$cov.mat <- cov.mat
results$std.errs <- std.errs
n.pars <- length(beta) + 2
results$bic <- (-2) * loglikhood + n.pars * log(sum(nvec))
results$aic <- (-2) * loglikhood + 2 * n.pars
return(results)
}
else {
N <- sum(nvec)
A1.st <- Diagonal(N, x = rep(1, N))
A2.st <- Diagonal(N, x = rep(1, N))
tmp <- BuildBasis(nvec)
E1.st <- tmp$E1
E2.st <- tmp$E2
E3.st <- tmp$E3
E4.st <- tmp$E4
R.alpha <- E1.st + E2.st + E3.st + E4.st
for (rr in 1:num.restarts) {
if (is.null(initvals) & rr == 1) {
initvals <- InitializePars(formula, data, subj.col,
nclasses, nvec, yy, XX, ylist, Xlist)
beta <- initvals$beta
alpha <- initvals$alpha
gamma <- initvals$gamma
mix.prop <- initvals$mix.prop
initvals.record[[1]] = initvals
}
else if (rr == 1 & !is.null(initvals)) {
beta <- initvals$coef
alpha <- initvals$autocorr
gamma <- initvals$scale - 1
mix.prop <- initvals$mix.prop
check.inits <- (class(beta) == "matrix") & (class(alpha) ==
"numeric") & (class(gamma) == "numeric") &
(class(mix.prop) == "numeric")
if (!check.inits) {
stop("Initial values are not in the correct format")
}
check.dims <- (nrow(beta) == nclasses) & (length(alpha) ==
nclasses) & (length(gamma) == nclasses) & (length(mix.prop) ==
nclasses)
if (!check.dims) {
stop("Initial values do not match the number of classes")
}
initvals.record[[1]] <- initvals
}
else if (rr > 1) {
invals <- MultiStart(parlist, nclasses, ylist,
Xlist)
beta <- invals$beta
alpha <- invals$alpha
gamma <- invals$gamma
mix.prop <- invals$mix.prop
initvals.record[[rr]] <- invals
}
newbeta <- matrix(0, nrow = nrow(beta), ncol = ncol(beta))
newalpha <- rep(0, length(alpha))
newgamma <- rep(0, length(gamma))
count <- 0
em.converged <- 0
current.max <- 1
GEE.conv <- matrix(0, nrow = maxiter + 1, ncol = nclasses)
loglikhood <- pss <- rep(0, maxiter + 1)
done <- FALSE
betalist <- lapply(1:nrow(beta), function(i) beta[i,
])
tmp.pp <- PostProbs(betalist, alpha, gamma, mix.prop,
ylist, Xlist, m, nclasses)
post.probs <- tmp.pp$postprobs
loglikhood[1] <- tmp.pp$loglik
for (i in 1:maxiter) {
if (maxiter == 0) {
count <- 1
break
}
for (j in 1:nclasses) {
initpars <- list(beta = beta[j, ], alpha = alpha[j],
gamma = gamma[j])
tmp <- GEEests(post.probs[j, ], initpars, XX,
yy, nvec, R.alpha, A1.st, A2.st, E1.st, E2.st,
E3.st, E4.st)
GEE.conv[i, j] <- tmp$conv
newparams <- EnforceConstraints(tmp, Xlist,
XX, subject.id)
newbeta[j, ] <- newparams$beta
newalpha[j] <- newparams$alpha
newgamma[j] <- newparams$gamma
}
newmix.prop <- rowSums(post.probs)/m
betalist <- lapply(1:nrow(newbeta), function(i) newbeta[i,
])
tmp.pp <- PostProbs(betalist, newalpha, newgamma,
newmix.prop, ylist, Xlist, m, nclasses)
post.probs <- tmp.pp$postprobs
loglikhood[i + 1] <- tmp.pp$loglik
par.stack <- c(c(newbeta), newalpha, newgamma,
newmix.prop[-nclasses])
pss.vec <- PsiFn(par.stack, len.beta = ncol(newbeta),
nclasses, yy, XX, ylist, Xlist, nvec, E1.st,
E2.st, E3.st, E4.st)
pss[i] <- sum(pss.vec^2)
if (i > 1) {
done <- all(abs(pss.vec) < stoptol)
}
if (is.na(done)) {
print("Warning: some post. probs are NA")
done <- FALSE
}
if (done) {
em.converged <- 1
count <- count + 1
beta <- newbeta
gamma <- newgamma
alpha <- newalpha
mix.prop <- newmix.prop
break
}
beta <- newbeta
gamma <- newgamma
alpha <- newalpha
mix.prop <- newmix.prop
if (any(mix.prop < dropthresh)) {
if (nclasses == 2) {
print("Try a one class model")
}
else {
print("Very small category, reducing the number of classes")
remove <- which(mix.prop == min(mix.prop))
mix.prop <- mix.prop[-remove]/sum(mix.prop[-remove])
nclasses <- nclasses - 1
beta <- as.matrix(beta[-remove, ])
alpha <- alpha[-remove]
gamma <- gamma[-remove]
newbeta <- matrix(0, nrow = nclasses, ncol = ncol(beta))
newalpha <- rep(0, nclasses)
newgamma <- rep(0, nclasses)
}
}
cat("Iteration: ", count + 1, "\n")
count <- count + 1
}
if (rr == 1) {
numiters <- count
EMconv <- em.converged
Loglik <- loglikhood[1:count]
GEEmat <- GEE.conv[1:count, ]
}
parlist <- list()
parlist$coef <- beta
parlist$auto.corr <- alpha
parlist$scale <- gamma + 1
parlist$mix.prop <- mix.prop
loglik.record[rr] <- loglikhood[count]
final.record[[rr]] <- parlist
if (rr > 1) {
if (loglik.record[rr] > loglik.record[rr - 1]) {
current.max <- rr
numiters <- count
EMconv <- em.converged
Loglik <- loglikhood[1:count]
GEEmat <- GEE.conv[1:count, ]
}
}
if (rr == num.restarts) {
beta <- final.record[[current.max]]$coef
alpha <- final.record[[current.max]]$auto.corr
gamma <- final.record[[current.max]]$scale -
1
mix.prop <- final.record[[current.max]]$mix.prop
reord <- order(-mix.prop)
beta <- as.matrix(beta[reord, ])
alpha <- alpha[reord]
gamma <- gamma[reord]
mix.prop <- mix.prop[reord]
betalist <- lapply(1:nrow(beta), function(i) beta[i,
])
tmp.pp <- PostProbs(betalist, alpha, gamma, mix.prop,
ylist, Xlist, m, nclasses)
post.probs <- tmp.pp$postprobs
q <- ncol(beta) + 2
par.stack <- numeric(nclasses * (q + 1) - 1)
for (i in 0:(nclasses - 1)) {
par.stack[(i * q + 1):((i + 1) * q)] <- c(beta[i +
1, ], alpha[i + 1], gamma[i + 1])
}
par.stack[(q * nclasses + 1):(nclasses * (q +
1) - 1)] <- mix.prop[-nclasses]
len.beta <- ncol(beta)
hh <- ComputeHess(par.stack, post.probs, len.beta,
ylist, Xlist, nvec)
if (rcond(hh$dermat) < sqrt(.Machine$double.eps)) {
print(rcond(hh$dermat))
hh$dermat <- hh$dermat + diag(rep(2e-08, nrow(hh$dermat)))
print(rcond(hh$dermat))
}
cov.mat <- tcrossprod(solve(hh$dermat, hh$sqmat),
solve(hh$dermat))
std.errs <- sqrt(diag(cov.mat))
tmp <- matrix("", ncol(XX) + 2, nclasses)
mix.names <- rep(c(""), nclasses)
for (k in 1:nclasses) {
for (j in 1:ncol(XX)) {
tmp[j, k] <- paste(colnames(XX)[j], k, sep = "")
}
tmp[ncol(XX) + 1, k] <- paste("autocorr.",
k, sep = "")
tmp[ncol(XX) + 2, k] <- paste("scale", k, sep = "")
mix.names[k] <- paste("mix", k)
}
names(std.errs) <- c(c(tmp), mix.names[-nclasses])
post.probs <- t(post.probs)
rownames(post.probs) <- subject.names
colnames(post.probs) <- rownames(beta) <- rep(" ",
nclasses)
colnames(beta) <- colnames(XX)
for (k in 1:nclasses) {
colnames(post.probs)[k] <- paste("group ",
k, sep = "")
rownames(beta)[k] <- paste("group ", k, sep = "")
}
results <- list()
class(results) <- "inarmix"
results$call <- mod
results$coefficients <- beta
results$alpha <- alpha
results$gamma <- gamma
results$mix.prop <- mix.prop
results$coefnames <- colnames(XX)
results$post.probs <- post.probs
results$fitted.values <- exp(tcrossprod(XX, beta))
res_se <- sqrt(t(t(results$fitted.values) * (gamma +
1)))
results$residuals <- (yy - results$fitted.values)/res_se
results$niter <- numiters
results$GEE.conv <- GEEmat
results$loglikhood <- Loglik
results$em.converged <- EMconv
results$pss <- pss[1:count]
results$initvals <- initvals.record[[current.max]]
results$nclasses <- nclasses
results$cov.mat <- cov.mat
results$std.errs <- std.errs
n.pars <- nclasses * (ncol(beta) + 3) - 1
results$bic <- (-2) * Loglik[numiters] + n.pars *
log(sum(nvec))
results$aic <- (-2) * Loglik[numiters] + 2 *
n.pars
if (num.restarts == 1) {
results$startingvals <- NULL
results$reploglik <- NULL
results$finalvals <- NULL
}
else {
results$startingvals <- initvals.record
results$reploglik <- loglik.record
results$finalvals <- final.record
}
}
else {
cat("\n Replication:", rr + 1, "\n")
}
}
return(results)
}
}
<bytecode: 0x63151c8>
<environment: namespace:inarmix>
--- function search by body ---
Function inarmix in namespace inarmix has this body.
----------- END OF FAILURE REPORT --------------
Error in if (!check.inits) { : the condition has length > 1
Calls: inarmix
Execution halted
Flavor: r-devel-linux-x86_64-debian-clang
Version: 0.4
Check: examples
Result: ERROR
Running examples in ‘inarmix-Ex.R’ failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: inarmix
> ### Title: Finite mixture model for longitudinal count data.
> ### Aliases: inarmix
> ### Keywords: models regression
>
> ### ** Examples
>
> set.seed(4297)
>
> ############################################################
> #### Simulate data from a two class model
>
> XX <- cbind(rep(1,9),c(0:8)/4)
> colnames(XX) <- c("const","time")
> beta <- rbind(c(-.2,0),c(1.2,.3))
> ### this means that for group 1: (beta_{0},beta_{1}) = (-.2,0)
> ### and for group 2: (beta_{0},beta_{1}) = (1.2,.3)
> autocorr <- c(.2,.2)
> scale <- c(2,2)
> mix.prop <- c(.8,.2) ## proportion in group 1 is .8
>
> testdat <- GenerateMixData(500,beta,autocorr,scale,mix.prop,XX)
> testdat[1:5,]
y subject const time
1 1 1 1 0.00
2 0 1 1 0.25
3 2 1 1 0.50
4 4 1 1 0.75
5 0 1 1 1.00
>
>
> ########################################################################
> #### Fit a linear curve with two classes (with a maximum of 4 iterations)
>
> twoclassfit <- inarmix(y~time,nclasses=2,id=subject,data=testdat,maxiter=4)
Iteration: 1
Iteration: 2
Iteration: 3
Iteration: 4
> summary(twoclassfit)
Call:
inarmix(formula = y ~ time, nclasses = 2, id = subject, data = testdat,
maxiter = 4)
Class Probabilities
Estimate Std.Err
Group1 0.7977567 0.01845
Group2 0.2022433 NA
------------------------------
Group 1 coefficients:
Estimate Std.Err Z value Pr(>z)
(Intercept) -0.18180172 0.05627233 -3.23075 0.0012347 **
time -0.02653325 0.04525100 -0.58636 0.5576354
autocorr. 0.20233773 0.02375121 8.51905 < 2.22e-16 ***
scale 2.00629134 0.08157575 12.33567 < 2.22e-16 ***
------------------------------
Group 2 coefficients:
Estimate Std.Err Z value Pr(>z)
(Intercept) 1.21738430 0.05813599 20.94029 < 2.22e-16 ***
time 0.28683640 0.04408047 6.50711 7.6611e-11 ***
autocorr. 0.19696018 0.03570303 5.51662 3.4557e-08 ***
scale 2.12021748 0.11034094 10.15233 < 2.22e-16 ***
------------------------------
log-likelihood: -6823.222
BIC: 13722.151
AIC: 13664.444
Number of EM iterations: 4
>
> diagnose(twoclassfit)
Did not converge!!!!
Convergence of GEE equations for
each iteration: 1 - converged, 0 - did not converge
iteration GEE 1 GEE 2 l[i] l[i+1]-l[i] ||Psi||^2
[1,] 0 1 1 -6866.134 0.000000 0.646708992
[2,] 1 1 1 -6833.577 -32.556835 0.108210548
[3,] 2 1 1 -6824.850 -8.726774 0.015731280
[4,] 3 1 1 -6823.222 -1.627973 0.002018329
>
> #############################################################
> ##### Fit the same model with specified starting values.
>
> inpars <- list()
> inpars$coef <- rbind(c(-.5,.1),c(.5,0))
> inpars$autocorr <- rep(.3,2)
> inpars$scale <- rep(2,2)
> inpars$mix.prop <- c(.6,.4)
>
> twoclassfit2 <- inarmix(y~time,nclasses=2,id=subject,data=testdat,initvals=inpars,
+ maxiter=4)
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
inarmix
--- call from context ---
inarmix(y ~ time, nclasses = 2, id = subject, data = testdat,
initvals = inpars, maxiter = 4)
--- call from argument ---
if (!check.inits) {
stop("Initial values are not in the correct format")
}
--- R stacktrace ---
where 1: inarmix(y ~ time, nclasses = 2, id = subject, data = testdat,
initvals = inpars, maxiter = 4)
--- value of length: 2 type: logical ---
[1] FALSE TRUE
--- function from context ---
function (formula, nclasses = 1, id, data, initvals = NULL, maxiter = 200,
stoptol = 1e-05, num.restarts = 1, time = NULL, dropthresh = 0.01)
{
mod <- match.call(expand.dots = FALSE)
modterms <- terms(formula)
subj.col <- which(colnames(data) == mod$id)
if (!is.null(mod$time)) {
time.col <- which(colnames(data) == mod$time)
if (length(time.col) == 0) {
stop("time variable not found")
}
}
if (length(subj.col) == 0) {
stop("id variable not found")
}
if (any(is.na(data))) {
stop("There are missing data points")
}
if (is.null(mod$time)) {
data <- data[order(data[, subj.col]), ]
}
else {
data <- data[order(data[, subj.col], data[, time.col]),
]
}
modframe <- model.frame(formula, data)
yy <- model.response(modframe)
XX <- model.matrix(formula, modframe)
QR <- qr(XX)
if (QR$rank < ncol(XX)) {
stop("The design matrix is rank-deficient")
}
if (any(is.infinite(yy)) | any(is.infinite(XX))) {
}
subject.id <- data[, subj.col]
m <- length(unique(subject.id))
nvec <- as.vector(table(subject.id))
loglik.record <- numeric(num.restarts)
initvals.record <- final.record <- list()
ylist <- Xlist <- list()
subject.names <- rep(" ", m)
ncum <- 1 + c(0, cumsum(nvec))
for (i in 1:m) {
ind1 <- ncum[i]
ind2 <- ncum[i + 1] - 1
ylist[[i]] <- yy[ind1:ind2]
Xlist[[i]] <- as.matrix(XX[ind1:ind2, ])
tmp <- unique(data[ind1, subj.col])
subject.names[i] <- as.character(tmp)
}
if (nclasses == 1) {
oneclass.ests <- OneclassEsts(yy, XX, nvec, m)
covs <- OuterProd(beta = oneclass.ests$beta, alpha = oneclass.ests$alpha,
gamma = oneclass.ests$gamma, nvec = nvec, ylist = ylist,
Xlist = Xlist)
cov.mat <- tcrossprod(solve(covs$D, covs$V), solve(covs$D))
std.errs <- sqrt(diag(cov.mat))
names(std.errs) <- c(colnames(XX), "autocorr.", "scale")
betalist <- list()
betalist[[1]] <- oneclass.ests$beta
tmp.pp <- PostProbs(betalist, oneclass.ests$alpha, oneclass.ests$gamma,
1, ylist, Xlist, m, 1)
post.probs <- tmp.pp[[1]]
loglikhood <- tmp.pp[[2]]
results <- list()
class(results) <- "inarmix"
results$call <- mod
results$coefficients <- oneclass.ests$beta
results$alpha <- oneclass.ests$alpha
results$gamma <- oneclass.ests$gamma
results$mix.prop <- NULL
results$coefnames <- colnames(XX)
results$nclasses <- 1
results$post.probs <- NULL
results$niter <- NULL
results$GEE.conv <- NULL
results$loglikhood <- loglikhood
results$em.converged <- NULL
results$pss <- NULL
results$initvals <- NULL
results$cov.mat <- cov.mat
results$std.errs <- std.errs
n.pars <- length(beta) + 2
results$bic <- (-2) * loglikhood + n.pars * log(sum(nvec))
results$aic <- (-2) * loglikhood + 2 * n.pars
return(results)
}
else {
N <- sum(nvec)
A1.st <- Diagonal(N, x = rep(1, N))
A2.st <- Diagonal(N, x = rep(1, N))
tmp <- BuildBasis(nvec)
E1.st <- tmp$E1
E2.st <- tmp$E2
E3.st <- tmp$E3
E4.st <- tmp$E4
R.alpha <- E1.st + E2.st + E3.st + E4.st
for (rr in 1:num.restarts) {
if (is.null(initvals) & rr == 1) {
initvals <- InitializePars(formula, data, subj.col,
nclasses, nvec, yy, XX, ylist, Xlist)
beta <- initvals$beta
alpha <- initvals$alpha
gamma <- initvals$gamma
mix.prop <- initvals$mix.prop
initvals.record[[1]] = initvals
}
else if (rr == 1 & !is.null(initvals)) {
beta <- initvals$coef
alpha <- initvals$autocorr
gamma <- initvals$scale - 1
mix.prop <- initvals$mix.prop
check.inits <- (class(beta) == "matrix") & (class(alpha) ==
"numeric") & (class(gamma) == "numeric") &
(class(mix.prop) == "numeric")
if (!check.inits) {
stop("Initial values are not in the correct format")
}
check.dims <- (nrow(beta) == nclasses) & (length(alpha) ==
nclasses) & (length(gamma) == nclasses) & (length(mix.prop) ==
nclasses)
if (!check.dims) {
stop("Initial values do not match the number of classes")
}
initvals.record[[1]] <- initvals
}
else if (rr > 1) {
invals <- MultiStart(parlist, nclasses, ylist,
Xlist)
beta <- invals$beta
alpha <- invals$alpha
gamma <- invals$gamma
mix.prop <- invals$mix.prop
initvals.record[[rr]] <- invals
}
newbeta <- matrix(0, nrow = nrow(beta), ncol = ncol(beta))
newalpha <- rep(0, length(alpha))
newgamma <- rep(0, length(gamma))
count <- 0
em.converged <- 0
current.max <- 1
GEE.conv <- matrix(0, nrow = maxiter + 1, ncol = nclasses)
loglikhood <- pss <- rep(0, maxiter + 1)
done <- FALSE
betalist <- lapply(1:nrow(beta), function(i) beta[i,
])
tmp.pp <- PostProbs(betalist, alpha, gamma, mix.prop,
ylist, Xlist, m, nclasses)
post.probs <- tmp.pp$postprobs
loglikhood[1] <- tmp.pp$loglik
for (i in 1:maxiter) {
if (maxiter == 0) {
count <- 1
break
}
for (j in 1:nclasses) {
initpars <- list(beta = beta[j, ], alpha = alpha[j],
gamma = gamma[j])
tmp <- GEEests(post.probs[j, ], initpars, XX,
yy, nvec, R.alpha, A1.st, A2.st, E1.st, E2.st,
E3.st, E4.st)
GEE.conv[i, j] <- tmp$conv
newparams <- EnforceConstraints(tmp, Xlist,
XX, subject.id)
newbeta[j, ] <- newparams$beta
newalpha[j] <- newparams$alpha
newgamma[j] <- newparams$gamma
}
newmix.prop <- rowSums(post.probs)/m
betalist <- lapply(1:nrow(newbeta), function(i) newbeta[i,
])
tmp.pp <- PostProbs(betalist, newalpha, newgamma,
newmix.prop, ylist, Xlist, m, nclasses)
post.probs <- tmp.pp$postprobs
loglikhood[i + 1] <- tmp.pp$loglik
par.stack <- c(c(newbeta), newalpha, newgamma,
newmix.prop[-nclasses])
pss.vec <- PsiFn(par.stack, len.beta = ncol(newbeta),
nclasses, yy, XX, ylist, Xlist, nvec, E1.st,
E2.st, E3.st, E4.st)
pss[i] <- sum(pss.vec^2)
if (i > 1) {
done <- all(abs(pss.vec) < stoptol)
}
if (is.na(done)) {
print("Warning: some post. probs are NA")
done <- FALSE
}
if (done) {
em.converged <- 1
count <- count + 1
beta <- newbeta
gamma <- newgamma
alpha <- newalpha
mix.prop <- newmix.prop
break
}
beta <- newbeta
gamma <- newgamma
alpha <- newalpha
mix.prop <- newmix.prop
if (any(mix.prop < dropthresh)) {
if (nclasses == 2) {
print("Try a one class model")
}
else {
print("Very small category, reducing the number of classes")
remove <- which(mix.prop == min(mix.prop))
mix.prop <- mix.prop[-remove]/sum(mix.prop[-remove])
nclasses <- nclasses - 1
beta <- as.matrix(beta[-remove, ])
alpha <- alpha[-remove]
gamma <- gamma[-remove]
newbeta <- matrix(0, nrow = nclasses, ncol = ncol(beta))
newalpha <- rep(0, nclasses)
newgamma <- rep(0, nclasses)
}
}
cat("Iteration: ", count + 1, "\n")
count <- count + 1
}
if (rr == 1) {
numiters <- count
EMconv <- em.converged
Loglik <- loglikhood[1:count]
GEEmat <- GEE.conv[1:count, ]
}
parlist <- list()
parlist$coef <- beta
parlist$auto.corr <- alpha
parlist$scale <- gamma + 1
parlist$mix.prop <- mix.prop
loglik.record[rr] <- loglikhood[count]
final.record[[rr]] <- parlist
if (rr > 1) {
if (loglik.record[rr] > loglik.record[rr - 1]) {
current.max <- rr
numiters <- count
EMconv <- em.converged
Loglik <- loglikhood[1:count]
GEEmat <- GEE.conv[1:count, ]
}
}
if (rr == num.restarts) {
beta <- final.record[[current.max]]$coef
alpha <- final.record[[current.max]]$auto.corr
gamma <- final.record[[current.max]]$scale -
1
mix.prop <- final.record[[current.max]]$mix.prop
reord <- order(-mix.prop)
beta <- as.matrix(beta[reord, ])
alpha <- alpha[reord]
gamma <- gamma[reord]
mix.prop <- mix.prop[reord]
betalist <- lapply(1:nrow(beta), function(i) beta[i,
])
tmp.pp <- PostProbs(betalist, alpha, gamma, mix.prop,
ylist, Xlist, m, nclasses)
post.probs <- tmp.pp$postprobs
q <- ncol(beta) + 2
par.stack <- numeric(nclasses * (q + 1) - 1)
for (i in 0:(nclasses - 1)) {
par.stack[(i * q + 1):((i + 1) * q)] <- c(beta[i +
1, ], alpha[i + 1], gamma[i + 1])
}
par.stack[(q * nclasses + 1):(nclasses * (q +
1) - 1)] <- mix.prop[-nclasses]
len.beta <- ncol(beta)
hh <- ComputeHess(par.stack, post.probs, len.beta,
ylist, Xlist, nvec)
if (rcond(hh$dermat) < sqrt(.Machine$double.eps)) {
print(rcond(hh$dermat))
hh$dermat <- hh$dermat + diag(rep(2e-08, nrow(hh$dermat)))
print(rcond(hh$dermat))
}
cov.mat <- tcrossprod(solve(hh$dermat, hh$sqmat),
solve(hh$dermat))
std.errs <- sqrt(diag(cov.mat))
tmp <- matrix("", ncol(XX) + 2, nclasses)
mix.names <- rep(c(""), nclasses)
for (k in 1:nclasses) {
for (j in 1:ncol(XX)) {
tmp[j, k] <- paste(colnames(XX)[j], k, sep = "")
}
tmp[ncol(XX) + 1, k] <- paste("autocorr.",
k, sep = "")
tmp[ncol(XX) + 2, k] <- paste("scale", k, sep = "")
mix.names[k] <- paste("mix", k)
}
names(std.errs) <- c(c(tmp), mix.names[-nclasses])
post.probs <- t(post.probs)
rownames(post.probs) <- subject.names
colnames(post.probs) <- rownames(beta) <- rep(" ",
nclasses)
colnames(beta) <- colnames(XX)
for (k in 1:nclasses) {
colnames(post.probs)[k] <- paste("group ",
k, sep = "")
rownames(beta)[k] <- paste("group ", k, sep = "")
}
results <- list()
class(results) <- "inarmix"
results$call <- mod
results$coefficients <- beta
results$alpha <- alpha
results$gamma <- gamma
results$mix.prop <- mix.prop
results$coefnames <- colnames(XX)
results$post.probs <- post.probs
results$fitted.values <- exp(tcrossprod(XX, beta))
res_se <- sqrt(t(t(results$fitted.values) * (gamma +
1)))
results$residuals <- (yy - results$fitted.values)/res_se
results$niter <- numiters
results$GEE.conv <- GEEmat
results$loglikhood <- Loglik
results$em.converged <- EMconv
results$pss <- pss[1:count]
results$initvals <- initvals.record[[current.max]]
results$nclasses <- nclasses
results$cov.mat <- cov.mat
results$std.errs <- std.errs
n.pars <- nclasses * (ncol(beta) + 3) - 1
results$bic <- (-2) * Loglik[numiters] + n.pars *
log(sum(nvec))
results$aic <- (-2) * Loglik[numiters] + 2 *
n.pars
if (num.restarts == 1) {
results$startingvals <- NULL
results$reploglik <- NULL
results$finalvals <- NULL
}
else {
results$startingvals <- initvals.record
results$reploglik <- loglik.record
results$finalvals <- final.record
}
}
else {
cat("\n Replication:", rr + 1, "\n")
}
}
return(results)
}
}
<bytecode: 0x5584ad4c1fd8>
<environment: namespace:inarmix>
--- function search by body ---
Function inarmix in namespace inarmix has this body.
----------- END OF FAILURE REPORT --------------
Error in if (!check.inits) { : the condition has length > 1
Calls: inarmix
Execution halted
Flavor: r-devel-linux-x86_64-debian-gcc
Version: 0.4
Check: compiled code
Result: NOTE
File ‘inarmix/libs/inarmix.so’:
Found no calls to: ‘R_registerRoutines’, ‘R_useDynamicSymbols’
It is good practice to register native routines and to disable symbol
search.
See ‘Writing portable packages’ in the ‘Writing R Extensions’ manual.
Flavors: r-devel-linux-x86_64-fedora-clang, r-devel-linux-x86_64-fedora-gcc
Version: 0.4
Check: examples
Result: ERROR
Running examples in ‘inarmix-Ex.R’ failed
The error most likely occurred in:
> ### Name: inarmix
> ### Title: Finite mixture model for longitudinal count data.
> ### Aliases: inarmix
> ### Keywords: models regression
>
> ### ** Examples
>
> set.seed(4297)
>
> ############################################################
> #### Simulate data from a two class model
>
> XX <- cbind(rep(1,9),c(0:8)/4)
> colnames(XX) <- c("const","time")
> beta <- rbind(c(-.2,0),c(1.2,.3))
> ### this means that for group 1: (beta_{0},beta_{1}) = (-.2,0)
> ### and for group 2: (beta_{0},beta_{1}) = (1.2,.3)
> autocorr <- c(.2,.2)
> scale <- c(2,2)
> mix.prop <- c(.8,.2) ## proportion in group 1 is .8
>
> testdat <- GenerateMixData(500,beta,autocorr,scale,mix.prop,XX)
> testdat[1:5,]
y subject const time
1 1 1 1 0.00
2 0 1 1 0.25
3 2 1 1 0.50
4 4 1 1 0.75
5 0 1 1 1.00
>
>
> ########################################################################
> #### Fit a linear curve with two classes (with a maximum of 4 iterations)
>
> twoclassfit <- inarmix(y~time,nclasses=2,id=subject,data=testdat,maxiter=4)
Iteration: 1
Iteration: 2
Iteration: 3
Iteration: 4
> summary(twoclassfit)
Call:
inarmix(formula = y ~ time, nclasses = 2, id = subject, data = testdat,
maxiter = 4)
Class Probabilities
Estimate Std.Err
Group1 0.7977567 0.01844
Group2 0.2022433 NA
------------------------------
Group 1 coefficients:
Estimate Std.Err Z value Pr(>z)
(Intercept) -0.18180172 0.05623476 -3.23291 0.0012254 **
time -0.02653325 0.04532190 -0.58544 0.5582519
autocorr. 0.20233773 0.02368218 8.54388 < 2.22e-16 ***
scale 2.00629134 0.08277197 12.15739 < 2.22e-16 ***
------------------------------
Group 2 coefficients:
Estimate Std.Err Z value Pr(>z)
(Intercept) 1.21738430 0.05829181 20.88431 < 2.22e-16 ***
time 0.28683640 0.04407928 6.50728 7.6522e-11 ***
autocorr. 0.19696018 0.03623996 5.43489 5.4830e-08 ***
scale 2.12021748 0.11447226 9.78593 < 2.22e-16 ***
------------------------------
log-likelihood: -6823.222
BIC: 13722.151
AIC: 13664.444
Number of EM iterations: 4
>
> diagnose(twoclassfit)
Did not converge!!!!
Convergence of GEE equations for
each iteration: 1 - converged, 0 - did not converge
iteration GEE 1 GEE 2 l[i] l[i+1]-l[i] ||Psi||^2
[1,] 0 1 1 -6866.134 0.000000 0.646708992
[2,] 1 1 1 -6833.577 -32.556835 0.108210548
[3,] 2 1 1 -6824.850 -8.726774 0.015731280
[4,] 3 1 1 -6823.222 -1.627973 0.002018329
>
> #############################################################
> ##### Fit the same model with specified starting values.
>
> inpars <- list()
> inpars$coef <- rbind(c(-.5,.1),c(.5,0))
> inpars$autocorr <- rep(.3,2)
> inpars$scale <- rep(2,2)
> inpars$mix.prop <- c(.6,.4)
>
> twoclassfit2 <- inarmix(y~time,nclasses=2,id=subject,data=testdat,initvals=inpars,
+ maxiter=4)
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
inarmix
--- call from context ---
inarmix(y ~ time, nclasses = 2, id = subject, data = testdat,
initvals = inpars, maxiter = 4)
--- call from argument ---
if (!check.inits) {
stop("Initial values are not in the correct format")
}
--- R stacktrace ---
where 1: inarmix(y ~ time, nclasses = 2, id = subject, data = testdat,
initvals = inpars, maxiter = 4)
--- value of length: 2 type: logical ---
[1] FALSE TRUE
--- function from context ---
function (formula, nclasses = 1, id, data, initvals = NULL, maxiter = 200,
stoptol = 1e-05, num.restarts = 1, time = NULL, dropthresh = 0.01)
{
mod <- match.call(expand.dots = FALSE)
modterms <- terms(formula)
subj.col <- which(colnames(data) == mod$id)
if (!is.null(mod$time)) {
time.col <- which(colnames(data) == mod$time)
if (length(time.col) == 0) {
stop("time variable not found")
}
}
if (length(subj.col) == 0) {
stop("id variable not found")
}
if (any(is.na(data))) {
stop("There are missing data points")
}
if (is.null(mod$time)) {
data <- data[order(data[, subj.col]), ]
}
else {
data <- data[order(data[, subj.col], data[, time.col]),
]
}
modframe <- model.frame(formula, data)
yy <- model.response(modframe)
XX <- model.matrix(formula, modframe)
QR <- qr(XX)
if (QR$rank < ncol(XX)) {
stop("The design matrix is rank-deficient")
}
if (any(is.infinite(yy)) | any(is.infinite(XX))) {
}
subject.id <- data[, subj.col]
m <- length(unique(subject.id))
nvec <- as.vector(table(subject.id))
loglik.record <- numeric(num.restarts)
initvals.record <- final.record <- list()
ylist <- Xlist <- list()
subject.names <- rep(" ", m)
ncum <- 1 + c(0, cumsum(nvec))
for (i in 1:m) {
ind1 <- ncum[i]
ind2 <- ncum[i + 1] - 1
ylist[[i]] <- yy[ind1:ind2]
Xlist[[i]] <- as.matrix(XX[ind1:ind2, ])
tmp <- unique(data[ind1, subj.col])
subject.names[i] <- as.character(tmp)
}
if (nclasses == 1) {
oneclass.ests <- OneclassEsts(yy, XX, nvec, m)
covs <- OuterProd(beta = oneclass.ests$beta, alpha = oneclass.ests$alpha,
gamma = oneclass.ests$gamma, nvec = nvec, ylist = ylist,
Xlist = Xlist)
cov.mat <- tcrossprod(solve(covs$D, covs$V), solve(covs$D))
std.errs <- sqrt(diag(cov.mat))
names(std.errs) <- c(colnames(XX), "autocorr.", "scale")
betalist <- list()
betalist[[1]] <- oneclass.ests$beta
tmp.pp <- PostProbs(betalist, oneclass.ests$alpha, oneclass.ests$gamma,
1, ylist, Xlist, m, 1)
post.probs <- tmp.pp[[1]]
loglikhood <- tmp.pp[[2]]
results <- list()
class(results) <- "inarmix"
results$call <- mod
results$coefficients <- oneclass.ests$beta
results$alpha <- oneclass.ests$alpha
results$gamma <- oneclass.ests$gamma
results$mix.prop <- NULL
results$coefnames <- colnames(XX)
results$nclasses <- 1
results$post.probs <- NULL
results$niter <- NULL
results$GEE.conv <- NULL
results$loglikhood <- loglikhood
results$em.converged <- NULL
results$pss <- NULL
results$initvals <- NULL
results$cov.mat <- cov.mat
results$std.errs <- std.errs
n.pars <- length(beta) + 2
results$bic <- (-2) * loglikhood + n.pars * log(sum(nvec))
results$aic <- (-2) * loglikhood + 2 * n.pars
return(results)
}
else {
N <- sum(nvec)
A1.st <- Diagonal(N, x = rep(1, N))
A2.st <- Diagonal(N, x = rep(1, N))
tmp <- BuildBasis(nvec)
E1.st <- tmp$E1
E2.st <- tmp$E2
E3.st <- tmp$E3
E4.st <- tmp$E4
R.alpha <- E1.st + E2.st + E3.st + E4.st
for (rr in 1:num.restarts) {
if (is.null(initvals) & rr == 1) {
initvals <- InitializePars(formula, data, subj.col,
nclasses, nvec, yy, XX, ylist, Xlist)
beta <- initvals$beta
alpha <- initvals$alpha
gamma <- initvals$gamma
mix.prop <- initvals$mix.prop
initvals.record[[1]] = initvals
}
else if (rr == 1 & !is.null(initvals)) {
beta <- initvals$coef
alpha <- initvals$autocorr
gamma <- initvals$scale - 1
mix.prop <- initvals$mix.prop
check.inits <- (class(beta) == "matrix") & (class(alpha) ==
"numeric") & (class(gamma) == "numeric") &
(class(mix.prop) == "numeric")
if (!check.inits) {
stop("Initial values are not in the correct format")
}
check.dims <- (nrow(beta) == nclasses) & (length(alpha) ==
nclasses) & (length(gamma) == nclasses) & (length(mix.prop) ==
nclasses)
if (!check.dims) {
stop("Initial values do not match the number of classes")
}
initvals.record[[1]] <- initvals
}
else if (rr > 1) {
invals <- MultiStart(parlist, nclasses, ylist,
Xlist)
beta <- invals$beta
alpha <- invals$alpha
gamma <- invals$gamma
mix.prop <- invals$mix.prop
initvals.record[[rr]] <- invals
}
newbeta <- matrix(0, nrow = nrow(beta), ncol = ncol(beta))
newalpha <- rep(0, length(alpha))
newgamma <- rep(0, length(gamma))
count <- 0
em.converged <- 0
current.max <- 1
GEE.conv <- matrix(0, nrow = maxiter + 1, ncol = nclasses)
loglikhood <- pss <- rep(0, maxiter + 1)
done <- FALSE
betalist <- lapply(1:nrow(beta), function(i) beta[i,
])
tmp.pp <- PostProbs(betalist, alpha, gamma, mix.prop,
ylist, Xlist, m, nclasses)
post.probs <- tmp.pp$postprobs
loglikhood[1] <- tmp.pp$loglik
for (i in 1:maxiter) {
if (maxiter == 0) {
count <- 1
break
}
for (j in 1:nclasses) {
initpars <- list(beta = beta[j, ], alpha = alpha[j],
gamma = gamma[j])
tmp <- GEEests(post.probs[j, ], initpars, XX,
yy, nvec, R.alpha, A1.st, A2.st, E1.st, E2.st,
E3.st, E4.st)
GEE.conv[i, j] <- tmp$conv
newparams <- EnforceConstraints(tmp, Xlist,
XX, subject.id)
newbeta[j, ] <- newparams$beta
newalpha[j] <- newparams$alpha
newgamma[j] <- newparams$gamma
}
newmix.prop <- rowSums(post.probs)/m
betalist <- lapply(1:nrow(newbeta), function(i) newbeta[i,
])
tmp.pp <- PostProbs(betalist, newalpha, newgamma,
newmix.prop, ylist, Xlist, m, nclasses)
post.probs <- tmp.pp$postprobs
loglikhood[i + 1] <- tmp.pp$loglik
par.stack <- c(c(newbeta), newalpha, newgamma,
newmix.prop[-nclasses])
pss.vec <- PsiFn(par.stack, len.beta = ncol(newbeta),
nclasses, yy, XX, ylist, Xlist, nvec, E1.st,
E2.st, E3.st, E4.st)
pss[i] <- sum(pss.vec^2)
if (i > 1) {
done <- all(abs(pss.vec) < stoptol)
}
if (is.na(done)) {
print("Warning: some post. probs are NA")
done <- FALSE
}
if (done) {
em.converged <- 1
count <- count + 1
beta <- newbeta
gamma <- newgamma
alpha <- newalpha
mix.prop <- newmix.prop
break
}
beta <- newbeta
gamma <- newgamma
alpha <- newalpha
mix.prop <- newmix.prop
if (any(mix.prop < dropthresh)) {
if (nclasses == 2) {
print("Try a one class model")
}
else {
print("Very small category, reducing the number of classes")
remove <- which(mix.prop == min(mix.prop))
mix.prop <- mix.prop[-remove]/sum(mix.prop[-remove])
nclasses <- nclasses - 1
beta <- as.matrix(beta[-remove, ])
alpha <- alpha[-remove]
gamma <- gamma[-remove]
newbeta <- matrix(0, nrow = nclasses, ncol = ncol(beta))
newalpha <- rep(0, nclasses)
newgamma <- rep(0, nclasses)
}
}
cat("Iteration: ", count + 1, "\n")
count <- count + 1
}
if (rr == 1) {
numiters <- count
EMconv <- em.converged
Loglik <- loglikhood[1:count]
GEEmat <- GEE.conv[1:count, ]
}
parlist <- list()
parlist$coef <- beta
parlist$auto.corr <- alpha
parlist$scale <- gamma + 1
parlist$mix.prop <- mix.prop
loglik.record[rr] <- loglikhood[count]
final.record[[rr]] <- parlist
if (rr > 1) {
if (loglik.record[rr] > loglik.record[rr - 1]) {
current.max <- rr
numiters <- count
EMconv <- em.converged
Loglik <- loglikhood[1:count]
GEEmat <- GEE.conv[1:count, ]
}
}
if (rr == num.restarts) {
beta <- final.record[[current.max]]$coef
alpha <- final.record[[current.max]]$auto.corr
gamma <- final.record[[current.max]]$scale -
1
mix.prop <- final.record[[current.max]]$mix.prop
reord <- order(-mix.prop)
beta <- as.matrix(beta[reord, ])
alpha <- alpha[reord]
gamma <- gamma[reord]
mix.prop <- mix.prop[reord]
betalist <- lapply(1:nrow(beta), function(i) beta[i,
])
tmp.pp <- PostProbs(betalist, alpha, gamma, mix.prop,
ylist, Xlist, m, nclasses)
post.probs <- tmp.pp$postprobs
q <- ncol(beta) + 2
par.stack <- numeric(nclasses * (q + 1) - 1)
for (i in 0:(nclasses - 1)) {
par.stack[(i * q + 1):((i + 1) * q)] <- c(beta[i +
1, ], alpha[i + 1], gamma[i + 1])
}
par.stack[(q * nclasses + 1):(nclasses * (q +
1) - 1)] <- mix.prop[-nclasses]
len.beta <- ncol(beta)
hh <- ComputeHess(par.stack, post.probs, len.beta,
ylist, Xlist, nvec)
if (rcond(hh$dermat) < sqrt(.Machine$double.eps)) {
print(rcond(hh$dermat))
hh$dermat <- hh$dermat + diag(rep(2e-08, nrow(hh$dermat)))
print(rcond(hh$dermat))
}
cov.mat <- tcrossprod(solve(hh$dermat, hh$sqmat),
solve(hh$dermat))
std.errs <- sqrt(diag(cov.mat))
tmp <- matrix("", ncol(XX) + 2, nclasses)
mix.names <- rep(c(""), nclasses)
for (k in 1:nclasses) {
for (j in 1:ncol(XX)) {
tmp[j, k] <- paste(colnames(XX)[j], k, sep = "")
}
tmp[ncol(XX) + 1, k] <- paste("autocorr.",
k, sep = "")
tmp[ncol(XX) + 2, k] <- paste("scale", k, sep = "")
mix.names[k] <- paste("mix", k)
}
names(std.errs) <- c(c(tmp), mix.names[-nclasses])
post.probs <- t(post.probs)
rownames(post.probs) <- subject.names
colnames(post.probs) <- rownames(beta) <- rep(" ",
nclasses)
colnames(beta) <- colnames(XX)
for (k in 1:nclasses) {
colnames(post.probs)[k] <- paste("group ",
k, sep = "")
rownames(beta)[k] <- paste("group ", k, sep = "")
}
results <- list()
class(results) <- "inarmix"
results$call <- mod
results$coefficients <- beta
results$alpha <- alpha
results$gamma <- gamma
results$mix.prop <- mix.prop
results$coefnames <- colnames(XX)
results$post.probs <- post.probs
results$fitted.values <- exp(tcrossprod(XX, beta))
res_se <- sqrt(t(t(results$fitted.values) * (gamma +
1)))
results$residuals <- (yy - results$fitted.values)/res_se
results$niter <- numiters
results$GEE.conv <- GEEmat
results$loglikhood <- Loglik
results$em.converged <- EMconv
results$pss <- pss[1:count]
results$initvals <- initvals.record[[current.max]]
results$nclasses <- nclasses
results$cov.mat <- cov.mat
results$std.errs <- std.errs
n.pars <- nclasses * (ncol(beta) + 3) - 1
results$bic <- (-2) * Loglik[numiters] + n.pars *
log(sum(nvec))
results$aic <- (-2) * Loglik[numiters] + 2 *
n.pars
if (num.restarts == 1) {
results$startingvals <- NULL
results$reploglik <- NULL
results$finalvals <- NULL
}
else {
results$startingvals <- initvals.record
results$reploglik <- loglik.record
results$finalvals <- final.record
}
}
else {
cat("\n Replication:", rr + 1, "\n")
}
}
return(results)
}
}
<bytecode: 0x5c089b8>
<environment: namespace:inarmix>
--- function search by body ---
Function inarmix in namespace inarmix has this body.
----------- END OF FAILURE REPORT --------------
Error in if (!check.inits) { : the condition has length > 1
Calls: inarmix
Execution halted
Flavor: r-devel-linux-x86_64-fedora-clang
Version: 0.4
Check: examples
Result: ERROR
Running examples in ‘inarmix-Ex.R’ failed
The error most likely occurred in:
> ### Name: inarmix
> ### Title: Finite mixture model for longitudinal count data.
> ### Aliases: inarmix
> ### Keywords: models regression
>
> ### ** Examples
>
> set.seed(4297)
>
> ############################################################
> #### Simulate data from a two class model
>
> XX <- cbind(rep(1,9),c(0:8)/4)
> colnames(XX) <- c("const","time")
> beta <- rbind(c(-.2,0),c(1.2,.3))
> ### this means that for group 1: (beta_{0},beta_{1}) = (-.2,0)
> ### and for group 2: (beta_{0},beta_{1}) = (1.2,.3)
> autocorr <- c(.2,.2)
> scale <- c(2,2)
> mix.prop <- c(.8,.2) ## proportion in group 1 is .8
>
> testdat <- GenerateMixData(500,beta,autocorr,scale,mix.prop,XX)
> testdat[1:5,]
y subject const time
1 1 1 1 0.00
2 0 1 1 0.25
3 2 1 1 0.50
4 4 1 1 0.75
5 0 1 1 1.00
>
>
> ########################################################################
> #### Fit a linear curve with two classes (with a maximum of 4 iterations)
>
> twoclassfit <- inarmix(y~time,nclasses=2,id=subject,data=testdat,maxiter=4)
Iteration: 1
Iteration: 2
Iteration: 3
Iteration: 4
> summary(twoclassfit)
Call:
inarmix(formula = y ~ time, nclasses = 2, id = subject, data = testdat,
maxiter = 4)
Class Probabilities
Estimate Std.Err
Group1 0.7977567 0.01845
Group2 0.2022433 NA
------------------------------
Group 1 coefficients:
Estimate Std.Err Z value Pr(>z)
(Intercept) -0.18180172 0.05627233 -3.23075 0.0012347 **
time -0.02653325 0.04525100 -0.58636 0.5576354
autocorr. 0.20233773 0.02375121 8.51905 < 2.22e-16 ***
scale 2.00629134 0.08157575 12.33567 < 2.22e-16 ***
------------------------------
Group 2 coefficients:
Estimate Std.Err Z value Pr(>z)
(Intercept) 1.21738430 0.05813599 20.94029 < 2.22e-16 ***
time 0.28683640 0.04408047 6.50711 7.6611e-11 ***
autocorr. 0.19696018 0.03570303 5.51662 3.4557e-08 ***
scale 2.12021748 0.11034094 10.15233 < 2.22e-16 ***
------------------------------
log-likelihood: -6823.222
BIC: 13722.151
AIC: 13664.444
Number of EM iterations: 4
>
> diagnose(twoclassfit)
Did not converge!!!!
Convergence of GEE equations for
each iteration: 1 - converged, 0 - did not converge
iteration GEE 1 GEE 2 l[i] l[i+1]-l[i] ||Psi||^2
[1,] 0 1 1 -6866.134 0.000000 0.646708992
[2,] 1 1 1 -6833.577 -32.556835 0.108210548
[3,] 2 1 1 -6824.850 -8.726774 0.015731280
[4,] 3 1 1 -6823.222 -1.627973 0.002018329
>
> #############################################################
> ##### Fit the same model with specified starting values.
>
> inpars <- list()
> inpars$coef <- rbind(c(-.5,.1),c(.5,0))
> inpars$autocorr <- rep(.3,2)
> inpars$scale <- rep(2,2)
> inpars$mix.prop <- c(.6,.4)
>
> twoclassfit2 <- inarmix(y~time,nclasses=2,id=subject,data=testdat,initvals=inpars,
+ maxiter=4)
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
inarmix
--- call from context ---
inarmix(y ~ time, nclasses = 2, id = subject, data = testdat,
initvals = inpars, maxiter = 4)
--- call from argument ---
if (!check.inits) {
stop("Initial values are not in the correct format")
}
--- R stacktrace ---
where 1: inarmix(y ~ time, nclasses = 2, id = subject, data = testdat,
initvals = inpars, maxiter = 4)
--- value of length: 2 type: logical ---
[1] FALSE TRUE
--- function from context ---
function (formula, nclasses = 1, id, data, initvals = NULL, maxiter = 200,
stoptol = 1e-05, num.restarts = 1, time = NULL, dropthresh = 0.01)
{
mod <- match.call(expand.dots = FALSE)
modterms <- terms(formula)
subj.col <- which(colnames(data) == mod$id)
if (!is.null(mod$time)) {
time.col <- which(colnames(data) == mod$time)
if (length(time.col) == 0) {
stop("time variable not found")
}
}
if (length(subj.col) == 0) {
stop("id variable not found")
}
if (any(is.na(data))) {
stop("There are missing data points")
}
if (is.null(mod$time)) {
data <- data[order(data[, subj.col]), ]
}
else {
data <- data[order(data[, subj.col], data[, time.col]),
]
}
modframe <- model.frame(formula, data)
yy <- model.response(modframe)
XX <- model.matrix(formula, modframe)
QR <- qr(XX)
if (QR$rank < ncol(XX)) {
stop("The design matrix is rank-deficient")
}
if (any(is.infinite(yy)) | any(is.infinite(XX))) {
}
subject.id <- data[, subj.col]
m <- length(unique(subject.id))
nvec <- as.vector(table(subject.id))
loglik.record <- numeric(num.restarts)
initvals.record <- final.record <- list()
ylist <- Xlist <- list()
subject.names <- rep(" ", m)
ncum <- 1 + c(0, cumsum(nvec))
for (i in 1:m) {
ind1 <- ncum[i]
ind2 <- ncum[i + 1] - 1
ylist[[i]] <- yy[ind1:ind2]
Xlist[[i]] <- as.matrix(XX[ind1:ind2, ])
tmp <- unique(data[ind1, subj.col])
subject.names[i] <- as.character(tmp)
}
if (nclasses == 1) {
oneclass.ests <- OneclassEsts(yy, XX, nvec, m)
covs <- OuterProd(beta = oneclass.ests$beta, alpha = oneclass.ests$alpha,
gamma = oneclass.ests$gamma, nvec = nvec, ylist = ylist,
Xlist = Xlist)
cov.mat <- tcrossprod(solve(covs$D, covs$V), solve(covs$D))
std.errs <- sqrt(diag(cov.mat))
names(std.errs) <- c(colnames(XX), "autocorr.", "scale")
betalist <- list()
betalist[[1]] <- oneclass.ests$beta
tmp.pp <- PostProbs(betalist, oneclass.ests$alpha, oneclass.ests$gamma,
1, ylist, Xlist, m, 1)
post.probs <- tmp.pp[[1]]
loglikhood <- tmp.pp[[2]]
results <- list()
class(results) <- "inarmix"
results$call <- mod
results$coefficients <- oneclass.ests$beta
results$alpha <- oneclass.ests$alpha
results$gamma <- oneclass.ests$gamma
results$mix.prop <- NULL
results$coefnames <- colnames(XX)
results$nclasses <- 1
results$post.probs <- NULL
results$niter <- NULL
results$GEE.conv <- NULL
results$loglikhood <- loglikhood
results$em.converged <- NULL
results$pss <- NULL
results$initvals <- NULL
results$cov.mat <- cov.mat
results$std.errs <- std.errs
n.pars <- length(beta) + 2
results$bic <- (-2) * loglikhood + n.pars * log(sum(nvec))
results$aic <- (-2) * loglikhood + 2 * n.pars
return(results)
}
else {
N <- sum(nvec)
A1.st <- Diagonal(N, x = rep(1, N))
A2.st <- Diagonal(N, x = rep(1, N))
tmp <- BuildBasis(nvec)
E1.st <- tmp$E1
E2.st <- tmp$E2
E3.st <- tmp$E3
E4.st <- tmp$E4
R.alpha <- E1.st + E2.st + E3.st + E4.st
for (rr in 1:num.restarts) {
if (is.null(initvals) & rr == 1) {
initvals <- InitializePars(formula, data, subj.col,
nclasses, nvec, yy, XX, ylist, Xlist)
beta <- initvals$beta
alpha <- initvals$alpha
gamma <- initvals$gamma
mix.prop <- initvals$mix.prop
initvals.record[[1]] = initvals
}
else if (rr == 1 & !is.null(initvals)) {
beta <- initvals$coef
alpha <- initvals$autocorr
gamma <- initvals$scale - 1
mix.prop <- initvals$mix.prop
check.inits <- (class(beta) == "matrix") & (class(alpha) ==
"numeric") & (class(gamma) == "numeric") &
(class(mix.prop) == "numeric")
if (!check.inits) {
stop("Initial values are not in the correct format")
}
check.dims <- (nrow(beta) == nclasses) & (length(alpha) ==
nclasses) & (length(gamma) == nclasses) & (length(mix.prop) ==
nclasses)
if (!check.dims) {
stop("Initial values do not match the number of classes")
}
initvals.record[[1]] <- initvals
}
else if (rr > 1) {
invals <- MultiStart(parlist, nclasses, ylist,
Xlist)
beta <- invals$beta
alpha <- invals$alpha
gamma <- invals$gamma
mix.prop <- invals$mix.prop
initvals.record[[rr]] <- invals
}
newbeta <- matrix(0, nrow = nrow(beta), ncol = ncol(beta))
newalpha <- rep(0, length(alpha))
newgamma <- rep(0, length(gamma))
count <- 0
em.converged <- 0
current.max <- 1
GEE.conv <- matrix(0, nrow = maxiter + 1, ncol = nclasses)
loglikhood <- pss <- rep(0, maxiter + 1)
done <- FALSE
betalist <- lapply(1:nrow(beta), function(i) beta[i,
])
tmp.pp <- PostProbs(betalist, alpha, gamma, mix.prop,
ylist, Xlist, m, nclasses)
post.probs <- tmp.pp$postprobs
loglikhood[1] <- tmp.pp$loglik
for (i in 1:maxiter) {
if (maxiter == 0) {
count <- 1
break
}
for (j in 1:nclasses) {
initpars <- list(beta = beta[j, ], alpha = alpha[j],
gamma = gamma[j])
tmp <- GEEests(post.probs[j, ], initpars, XX,
yy, nvec, R.alpha, A1.st, A2.st, E1.st, E2.st,
E3.st, E4.st)
GEE.conv[i, j] <- tmp$conv
newparams <- EnforceConstraints(tmp, Xlist,
XX, subject.id)
newbeta[j, ] <- newparams$beta
newalpha[j] <- newparams$alpha
newgamma[j] <- newparams$gamma
}
newmix.prop <- rowSums(post.probs)/m
betalist <- lapply(1:nrow(newbeta), function(i) newbeta[i,
])
tmp.pp <- PostProbs(betalist, newalpha, newgamma,
newmix.prop, ylist, Xlist, m, nclasses)
post.probs <- tmp.pp$postprobs
loglikhood[i + 1] <- tmp.pp$loglik
par.stack <- c(c(newbeta), newalpha, newgamma,
newmix.prop[-nclasses])
pss.vec <- PsiFn(par.stack, len.beta = ncol(newbeta),
nclasses, yy, XX, ylist, Xlist, nvec, E1.st,
E2.st, E3.st, E4.st)
pss[i] <- sum(pss.vec^2)
if (i > 1) {
done <- all(abs(pss.vec) < stoptol)
}
if (is.na(done)) {
print("Warning: some post. probs are NA")
done <- FALSE
}
if (done) {
em.converged <- 1
count <- count + 1
beta <- newbeta
gamma <- newgamma
alpha <- newalpha
mix.prop <- newmix.prop
break
}
beta <- newbeta
gamma <- newgamma
alpha <- newalpha
mix.prop <- newmix.prop
if (any(mix.prop < dropthresh)) {
if (nclasses == 2) {
print("Try a one class model")
}
else {
print("Very small category, reducing the number of classes")
remove <- which(mix.prop == min(mix.prop))
mix.prop <- mix.prop[-remove]/sum(mix.prop[-remove])
nclasses <- nclasses - 1
beta <- as.matrix(beta[-remove, ])
alpha <- alpha[-remove]
gamma <- gamma[-remove]
newbeta <- matrix(0, nrow = nclasses, ncol = ncol(beta))
newalpha <- rep(0, nclasses)
newgamma <- rep(0, nclasses)
}
}
cat("Iteration: ", count + 1, "\n")
count <- count + 1
}
if (rr == 1) {
numiters <- count
EMconv <- em.converged
Loglik <- loglikhood[1:count]
GEEmat <- GEE.conv[1:count, ]
}
parlist <- list()
parlist$coef <- beta
parlist$auto.corr <- alpha
parlist$scale <- gamma + 1
parlist$mix.prop <- mix.prop
loglik.record[rr] <- loglikhood[count]
final.record[[rr]] <- parlist
if (rr > 1) {
if (loglik.record[rr] > loglik.record[rr - 1]) {
current.max <- rr
numiters <- count
EMconv <- em.converged
Loglik <- loglikhood[1:count]
GEEmat <- GEE.conv[1:count, ]
}
}
if (rr == num.restarts) {
beta <- final.record[[current.max]]$coef
alpha <- final.record[[current.max]]$auto.corr
gamma <- final.record[[current.max]]$scale -
1
mix.prop <- final.record[[current.max]]$mix.prop
reord <- order(-mix.prop)
beta <- as.matrix(beta[reord, ])
alpha <- alpha[reord]
gamma <- gamma[reord]
mix.prop <- mix.prop[reord]
betalist <- lapply(1:nrow(beta), function(i) beta[i,
])
tmp.pp <- PostProbs(betalist, alpha, gamma, mix.prop,
ylist, Xlist, m, nclasses)
post.probs <- tmp.pp$postprobs
q <- ncol(beta) + 2
par.stack <- numeric(nclasses * (q + 1) - 1)
for (i in 0:(nclasses - 1)) {
par.stack[(i * q + 1):((i + 1) * q)] <- c(beta[i +
1, ], alpha[i + 1], gamma[i + 1])
}
par.stack[(q * nclasses + 1):(nclasses * (q +
1) - 1)] <- mix.prop[-nclasses]
len.beta <- ncol(beta)
hh <- ComputeHess(par.stack, post.probs, len.beta,
ylist, Xlist, nvec)
if (rcond(hh$dermat) < sqrt(.Machine$double.eps)) {
print(rcond(hh$dermat))
hh$dermat <- hh$dermat + diag(rep(2e-08, nrow(hh$dermat)))
print(rcond(hh$dermat))
}
cov.mat <- tcrossprod(solve(hh$dermat, hh$sqmat),
solve(hh$dermat))
std.errs <- sqrt(diag(cov.mat))
tmp <- matrix("", ncol(XX) + 2, nclasses)
mix.names <- rep(c(""), nclasses)
for (k in 1:nclasses) {
for (j in 1:ncol(XX)) {
tmp[j, k] <- paste(colnames(XX)[j], k, sep = "")
}
tmp[ncol(XX) + 1, k] <- paste("autocorr.",
k, sep = "")
tmp[ncol(XX) + 2, k] <- paste("scale", k, sep = "")
mix.names[k] <- paste("mix", k)
}
names(std.errs) <- c(c(tmp), mix.names[-nclasses])
post.probs <- t(post.probs)
rownames(post.probs) <- subject.names
colnames(post.probs) <- rownames(beta) <- rep(" ",
nclasses)
colnames(beta) <- colnames(XX)
for (k in 1:nclasses) {
colnames(post.probs)[k] <- paste("group ",
k, sep = "")
rownames(beta)[k] <- paste("group ", k, sep = "")
}
results <- list()
class(results) <- "inarmix"
results$call <- mod
results$coefficients <- beta
results$alpha <- alpha
results$gamma <- gamma
results$mix.prop <- mix.prop
results$coefnames <- colnames(XX)
results$post.probs <- post.probs
results$fitted.values <- exp(tcrossprod(XX, beta))
res_se <- sqrt(t(t(results$fitted.values) * (gamma +
1)))
results$residuals <- (yy - results$fitted.values)/res_se
results$niter <- numiters
results$GEE.conv <- GEEmat
results$loglikhood <- Loglik
results$em.converged <- EMconv
results$pss <- pss[1:count]
results$initvals <- initvals.record[[current.max]]
results$nclasses <- nclasses
results$cov.mat <- cov.mat
results$std.errs <- std.errs
n.pars <- nclasses * (ncol(beta) + 3) - 1
results$bic <- (-2) * Loglik[numiters] + n.pars *
log(sum(nvec))
results$aic <- (-2) * Loglik[numiters] + 2 *
n.pars
if (num.restarts == 1) {
results$startingvals <- NULL
results$reploglik <- NULL
results$finalvals <- NULL
}
else {
results$startingvals <- initvals.record
results$reploglik <- loglik.record
results$finalvals <- final.record
}
}
else {
cat("\n Replication:", rr + 1, "\n")
}
}
return(results)
}
}
<bytecode: 0x7709658>
<environment: namespace:inarmix>
--- function search by body ---
Function inarmix in namespace inarmix has this body.
----------- END OF FAILURE REPORT --------------
Error in if (!check.inits) { : the condition has length > 1
Calls: inarmix
Execution halted
Flavor: r-devel-linux-x86_64-fedora-gcc
Version: 0.4
Check: running examples for arch ‘i386’
Result: ERROR
Running examples in 'inarmix-Ex.R' failed
The error most likely occurred in:
> ### Name: diagnose
> ### Title: Diagnostics for the model fit.
> ### Aliases: diagnose
> ### Keywords: print
>
> ### ** Examples
>
> XX <- cbind(rep(1,9),c(0:8)/4)
> colnames(XX) <- c("const","time")
> coefs <- rbind(c(-.2,0),c(1.2,.3))
> alpha <- c(.2,.2)
> scale <- c(2,2)
> mix.prop <- c(.8,.2)
>
> testdat <- GenerateMixData(200,coefs,alpha,scale,mix.prop,XX)
> testfit <- inarmix(y~time,nclasses=2,id=subject,data=testdat,maxiter=3)
Iteration: 1
Iteration: 2
Iteration: 3
[1] 8.899245e-157
[1] 8.899245e-157
Error in solve.default(hh$dermat, hh$sqmat) :
system is computationally singular: reciprocal condition number = 8.89925e-157
Calls: inarmix -> tcrossprod -> solve -> solve -> solve.default
Execution halted
Flavor: r-devel-windows-ix86+x86_64-gcc8