File: JointAutoStrat.R

package info (click to toggle)
r-cran-openmx 2.21.1%2Bdfsg-1
  • links: PTS, VCS
  • area: main
  • in suites: bookworm
  • size: 14,412 kB
  • sloc: cpp: 36,577; ansic: 13,811; fortran: 2,001; sh: 1,440; python: 350; perl: 21; makefile: 5
file content (74 lines) | stat: -rw-r--r-- 2,347 bytes parent folder | download | duplicates (3)
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
libraries <- rownames(installed.packages())
if (!all(c("emx","gmp") %in% libraries)) stop("SKIP")

library(OpenMx)
library(emx)
library(rpf)
library(gmp)

numPatInclude <- c()
for (np in seq(2,200,2)) {
  numOutcome <- factorize(np)
  if (any(numOutcome > 7)) next
  numPatInclude <- c(numPatInclude, np)
}

todo <- expand.grid(numPat=numPatInclude, numContinuous=c(1L,100L),
                    numUnique=NA, strat=NA)

for (trial in 1:nrow(todo)) {
#for (trial in which(is.na(todo$continuous))) {
  set.seed(trial)
  numPat <- todo[trial,'numPat']
  exampleData <- list()
  numOutcome <- factorize(numPat)
  numOrdinal <- length(numOutcome)
  for (ix in 1:numOrdinal) {
    exampleData[[paste0('o',ix)]] <-
      mxFactor(1, levels=1:as.integer(numOutcome[ix]))
  }
  exampleData <- as.data.frame(exampleData)
  
  numContinuous <- todo[trial,'numContinuous']
  
  manifests <- c(paste0("o", 1:numOrdinal), paste0("c", 1:numContinuous))
  
  j1 <- mxModel(
    'j1', type="RAM",
    manifestVars=manifests,
    latentVars = 'G',
    emxThresholds(exampleData),
    mxData(exampleData, 'raw'),
    mxPath('G', manifests, values=runif(length(manifests), .5,1.5)),
    mxPath('G', free=FALSE, arrows=2, values=1),
    mxPath(paste0("o", 1:numOrdinal), arrows=2,
           values=.1+rlnorm(numOrdinal), free=FALSE),
    mxPath(paste0("c", 1:numContinuous), arrows=2,
           values=.1+rlnorm(numContinuous), free=TRUE),
    mxPath('one', paste0("c", 1:numContinuous),
           values=runif(numContinuous,-.5,.5)),
    mxComputeOnce('fitfunction','fit'))
  
  j1$expectation$thresholds <- "thresholdMatrix"
  
  ss <- 1000
  j2 <- mxGenerateData(j1, nrows=ss, returnModel=TRUE)
  cdf <- compressDataFrame(j2$data$observed[,paste0('o', 1:numOrdinal),drop=FALSE])
  todo[trial,'numUnique'] <- nrow(cdf)

  if (j2$fitfunction$jointConditionOn != 'auto') stop("not auto?")

#  j2$fitfunction$verbose <- 2L
  got <- mxRun(j2, silent = TRUE)

  strat <- attr(got$fitfunction$result,'jointConditionOn')
  if (!is.factor(todo$strat)) {
    todo$strat <- factor(todo$strat, levels=levels(strat))
  }
  todo[trial,c('strat')] <- strat
}

omxCheckTrue(todo$numUnique <= todo$numPat)
pred <- -3.58 + 0.36 * 1000/todo$numUnique - 0.06 * todo$numContinuous
omxCheckEquals(todo[pred<0, 'strat'], 'continuous')
omxCheckEquals(todo[pred>0, 'strat'], 'ordinal')