1char2num <- function(x,...) {
2    idx <- grep("^[-]*[0-9\\.]+",x,perl=TRUE,invert=TRUE)
3    if (length(idx)>0) x[idx] <- NA
4    as.numeric(x)
5}
6
7###{{{ substArg
8
9substArg <- function(x,env,...) {
10  if (!missing(env)) {
11    a <- with(env,substitute(x))
12  } else {
13    a <- substitute(x)
14  }
15  myclass <- tryCatch(class(eval(a)),error=function(e) NULL)
16  if (is.null(myclass) || myclass=="name") {
17      res <- unlist(sapply(as.character(a),
18                           function(z) {
19                               trimmed <- gsub(" ","",z,fixed=TRUE)
20                               val <- strsplit(trimmed,"+",fixed=TRUE)
21                               if (val[1]=="") val <- NULL
22                               val
23                           })); attributes(res)$names <- NULL
24      return(res)
25  }
26  return(eval(a))
27}
28
29## g <- function(zz,...) {
30##   env=new.env(); assign("x",substitute(zz),env)
31##   substArg(zz,env=env)
32## }
33## h <- function(x,...) {
34##   env=new.env(); assign("x",substitute(x),env)
35##   substArg(x,env=TRUE)
36## }
37
38###}}}
39
40###{{{ procrandomslope
41
42procrandomslope <- function(object,data=object$data,...) {
43  Xfix <- FALSE
44  xfix <- myfix <- list()
45  xx <- object
46  for (i in seq_len(object$ngroup)) {
47    x0 <- object$lvm[[i]]
48    data0 <- data[[i]]
49    xfix0 <- colnames(data0)[(colnames(data0)%in%parlabels(x0,exo=TRUE))]
50    xfix <- c(xfix, list(xfix0))
51    if (length(xfix0)>0) { ## Yes, random slopes
52      Xfix<-TRUE
53    }
54    xx$lvm[[i]] <- x0
55  }
56  if (Xfix) {
57    for (k in seq_len(object$ngroup)) {
58      x0 <- object$lvm[[k]]
59      data0 <- data[[k]]
60      nrow <- length(vars(x0))
61      xpos <- lapply(xfix[[k]],function(y) which(regfix(x0)$labels==y))
62      colpos <- lapply(xpos, function(y) ceiling(y/nrow))
63      rowpos <- lapply(xpos, function(y) (y-1)%%nrow+1)
64      myfix0 <- list(var=xfix[[k]], col=colpos, row=rowpos)
65      myfix <- c(myfix, list(myfix0))
66      for (i in seq_along(myfix0$var))
67        for (j in seq_along(myfix0$col[[i]]))
68          regfix(x0,
69                 from=vars(x0)[myfix0$row[[i]][j]],to=vars(x0)[myfix0$col[[i]][j]]) <-
70                   colMeans(data0[,myfix0$var[[i]],drop=FALSE],na.rm=TRUE)
71      index(x0) <- reindex(x0,zeroones=TRUE,deriv=TRUE)
72      object$lvm[[k]] <- x0
73    }
74    object <- multigroup(object$lvm,data,fix=FALSE,exo.fix=FALSE)
75  }
76  return(list(model=object,fix=myfix))
77}
78
79###}}} procrandomslope
80
81###{{{ kronprod
82
83## ' Calculate matrix product with kronecker product
84## '
85## ' \deqn{(A\crossprod B) Y}
86## ' @title Calculate matrix product with kronecker product
87## ' @param A
88## ' @param B
89## ' @param Y
90## ' @author Klaus K. Holst
91kronprod <- function(A,B,Y) {
92    if (missing(Y)) {
93        ## Assume 'B'=Identity, (A otimes B)Y
94        k <- nrow(B)/ncol(A)
95        res <- rbind(apply(B,2,function(x) matrix(x,nrow=k)%*%t(A)))
96        return(res)
97    }
98    rbind(apply(Y,2,function(x) B%*%matrix(x,nrow=ncol(B))%*%t(A)))
99}
100
101###}}} kronprod
102
103###{{{ izero
104
105izero <- function(i,n) { ## n-1 zeros and 1 at ith entry
106  x <- rep(0,n); x[i] <- 1
107  x
108}
109
110###}}}
111
112###{{{ Debug
113
114`Debug` <-
115  function(msg, cond=lava.options()$debug) {
116    if (cond)
117      print(paste(msg, collapse=" "))
118  }
119
120###}}}
121
122###{{{ categorical2dummy
123
124categorical2dummy <- function(x,data,messages=0,...) {
125  x0 <- x
126  X <- intersect(index(x)$exogenous,colnames(data))
127  catX <- c()
128  for (i in X) {
129    if (!is.numeric(data[,i])) catX <- c(catX,i)
130  }
131  if (length(catX)==0) return(list(x=x,data=data))
132  f <- as.formula(paste("~ 1+", paste(catX,collapse="+")))
133  opt <- options(na.action="na.pass")
134  M <- model.matrix(f,data)
135
136  options(opt)
137  Mnames <- colnames(M)
138  Mpos <- attributes(M)$assign
139  A <- index(x)$A
140  F <- regfix(x)
141  count <- 0
142  for (i in catX) {
143    count <- count+1
144    mnames <- Mnames[Mpos==count]
145    kill(x0) <- i
146    Y <- colnames(A)[A[i,]==1]
147    if (length(mnames)==1) {
148      fix <- as.list(F$labels[i,])
149      fixval <- F$values[i,]
150      fix[which(!is.na(fixval))] <- fixval[na.omit(fixval)]
151      regression(x0,to=Y,from=mnames,messages=messages) <- fix[Y]
152    } else {
153      x0 <- regression(x0,to=Y,from=mnames,messages=messages)
154    }
155  }
156  index(x0) <- reindex(x0,zeroones=TRUE,deriv=TRUE)
157  return(list(x=x0,data=cbind(data,M)))
158}
159
160###}}}
161
162###{{{ procdata.lvm
163
164`procdata.lvm` <-
165  function(x,data,categorical=FALSE,
166    na.method=ifelse(any(is.na(data[,intersect(colnames(data),manifest(x))])),"complete.obs","pairwise.complete.obs"),
167    missing=FALSE
168    ) {
169    if (is.numeric(data) & !is.list(data)) {
170      data <- rbind(data)
171    }
172     if (is.data.frame(data) | is.matrix(data)) {
173      nn <- colnames(data)
174      data <- as.data.frame(data); colnames(data) <- nn; rownames(data) <- NULL
175      obs <- setdiff(intersect(vars(x), colnames(data)),latent(x))
176      Debug(obs)
177      mydata <- subset(data, select=obs)
178      if (NROW(mydata)==0) stop("No observations")
179      for (i in seq_len(ncol(mydata))) {
180        if (inherits(mydata[,i],"Surv"))
181          mydata[,i] <- mydata[,i][,1]
182        if (is.character(mydata[,i]) | is.factor(mydata[,i]))
183          mydata[,i] <- as.numeric(as.factor(mydata[,i]))-1
184      }
185
186      S <- NULL
187      n <- nrow(mydata)
188      if (n==1) {
189        S <- diag(nrow=ncol(mydata)); colnames(S) <- rownames(S) <- obs
190      }
191      if (na.method=="complete.obs" && !missing) {
192        mydata0 <- na.omit(mydata)
193        n <- nrow(mydata0)
194        mu <- colMeans(mydata0)
195        if (is.null(S) && n>2)
196            S <- (n-1)/n*cov(mydata0) ## MLE variance matrix of observed variables
197        rm(mydata0)
198      }
199      nS <- is.null(S) || any(is.na(S))
200      if (na.method=="pairwise.complete.obs" || nS) {
201          mu <- colMeans(mydata,na.rm=TRUE)
202          if (nS) {
203              n <- nrow(mydata)
204              S <- (n-1)/n*cov(mydata,use="pairwise.complete.obs")
205              S[is.na(S)] <- 1e-3
206          }
207      }
208    }
209    else
210      if (is.list(data)) {
211        if ("cov"%in%names(data)) data$S <- data$cov
212        if ("var"%in%names(data)) data$S <- data$var
213        if ("mean"%in%names(data)) data$mu <- data$mean
214        n <- data$n
215        S <- reorderdata.lvm(x,data$S)
216        mu <- reorderdata.lvm(x,data$mu)
217        ##      if (is.null(n)) stop("n was not specified");
218      }
219      else
220        stop("Unexpected type of data!");
221    if (nrow(S)!=ncol(S)) stop("Wrong type of data!");
222    return(list(S=S,mu=mu,n=n))
223  }
224
225###}}}
226
227###{{{ reorderdata.lvm
228
229`reorderdata.lvm` <-
230  function(x, data) {
231    if (is.vector(data)) {
232      nn <- names(data)
233      ii <- na.omit(match(index(x)$manifest, nn))
234      data[ii,drop=FALSE]
235    } else {
236      nn <- colnames(data)
237      ii <- na.omit(match(index(x)$manifest, nn))
238      data[ii,ii,drop=FALSE]
239    }
240  }
241
242###}}}
243
244###{{{ symmetrize
245
246`symmetrize` <-
247function(M, upper=TRUE) {
248  if (length(M)==1) return(M)
249  if (!is.matrix(M) | ncol(M)!=nrow(M)) stop("Only implemented for square matrices.")
250  if (upper) {
251    for (i in seq_len(ncol(M)-1))
252      for (j in seq(i+1,nrow(M)))
253        M[i,j] <- M[j,i]
254    return(M)
255  } else {
256    for (i in seq_len(ncol(M)))
257      for (j in seq_len(nrow(M)))
258        if (M[i,j]==0)
259          M[i,j] <- M[j,i]
260        else
261          M[j,i] <- M[i,j]
262    return(M)
263  }
264}
265
266###}}}
267
268###{{{ naiveGrad
269
270naiveGrad <- function(f, x, h=1e-9) {
271  nabla <- numeric(length(x))
272  for (i in seq_along(x)) {
273    xh <- x; xh[i] <- x[i]+h
274    nabla[i] <- (f(xh)-f(x))/h
275  }
276  return(nabla)
277}
278
279###}}}
280
281###{{{ CondMom
282
283# conditional on Compl(idx)
284CondMom <- function(mu,S,idx,X) {
285  idxY <- idx
286
287  idxX <- setdiff(seq_len(ncol(S)),idxY)
288  SXX <- S[idxX,idxX,drop=FALSE];
289  SYY <- S[idxY,idxY,drop=FALSE]
290  SYX <- S[idxY,idxX,drop=FALSE]
291  iSXX <- solve(SXX)
292  condvar <- SYY-SYX%*%iSXX%*%t(SYX)
293  if (missing(mu)) return(condvar)
294
295  muY <- mu[,idxY,drop=FALSE]
296  muX <- mu[,idxX,drop=FALSE]
297  if (is.matrix(mu))
298    Z <- t(X-muX)
299  else
300    Z <- apply(X,1,function(xx) xx-muX)
301  SZ  <- t(SYX%*%iSXX%*%Z)
302##  condmean <- matrix(
303  if (is.matrix(mu))
304    condmean <- SZ+muY
305  else
306    condmean <- t(apply(SZ,1,function(x) muY+x))
307##  ,ncol=ncol(SZ),nrow=nrow(SZ))
308  return(list(mean=condmean,var=condvar))
309}
310
311###}}} CondMom
312
313###{{{ Depth-First/acc (accessible)
314
315DFS <- function(M,v,explored=c()) {
316  explored <- union(explored,v)
317  incident <- M[v,]
318  for (v1 in setdiff(which(incident==1),explored)) {
319    explored <- DFS(M,v1,explored)
320  }
321  return(explored)
322}
323
324acc <- function(M,v) {
325  if (is.character(v)) v <- which(colnames(M)==v)
326  colnames(M)[setdiff(DFS(M,v),v)]
327}
328
329###}}} Depth-First/acc (accessible)
330
331
332npar.lvm <- function(x) {
333  return(index(x)$npar+ index(x)$npar.mean+index(x)$npar.ex)
334
335}
336
337as.numeric.list <- function(x,...) {
338  lapply(x,function(y) ifelse(is.na(as.numeric(y)),y,as.numeric(y)))
339}
340
341edge2pair <- function(e) {
342  sapply(e,function(x) strsplit(x,"~"))
343}
344numberdup <- function(xx) { ## Convert to numbered list
345  dup.xx <- duplicated(xx)
346  ## dups <- xx[dup.xx]
347  xx.new <- numeric(length(xx))
348  count <- 0
349  for (i in seq_along(xx)) {
350    if (!dup.xx[i]) {
351      count <- count+1
352      xx.new[i] <- count
353    } else {
354      xx.new[i] <- xx.new[match(xx[i],xx)[1]]
355    }
356  }
357  return(xx.new)
358}
359
360extractvar <- function(f) {
361    yy <- getoutcome(f)
362    xx <- attributes(terms(f))$term.labels
363    myvars <- all.vars(f)
364    return(list(y=yy,x=xx,all=myvars))
365}
366
367##' @export
368getoutcome <- function(formula,sep,...) {
369  aa <- attributes(terms(formula,...))
370  if (aa$response==0) {
371    res <- NULL
372  } else {
373    res <- paste(deparse(formula[[2]]),collapse="")
374  }
375  if (!missing(sep) && length(aa$term.labels)>0) {
376      attributes(res)$x <- lapply(strsplit(aa$term.labels,"\\|")[[1]],
377                                  function(x) as.formula(paste0("~",x)))
378  } else {
379      attributes(res)$x <- aa$term.labels
380  }
381  return(res)
382}
383
384
385##' @export
386Specials <- function(f,spec,split2="+",...) {
387  tt <- terms(f,spec)
388  pos <- attributes(tt)$specials[[spec]]
389  if (is.null(pos)) return(NULL)
390  x <- rownames(attributes(tt)$factors)[pos]
391  st <- gsub(" ","",x)
392  res <- unlist(strsplit(st,"[()]"))[2]
393  if (is.null(split2)) return(res)
394  unlist(strsplit(res,split2,fixed=TRUE))
395}
396
397
398##' @export
399decomp.specials <- function(x,pattern="[()]",pattern2=NULL, pattern.ignore=NULL, sep="[,\\+]",perl=TRUE,reverse=FALSE,...) {
400  st <- gsub(" |^\\(|)$","",x) # Remove white space and leading/trailing parantheses
401  if (!is.null(pattern.ignore)) {
402      if (grepl(pattern.ignore,st,perl=perl,...)) return(st)
403  }
404  if (!is.null(pattern)) {
405    st <- rev(unlist(strsplit(st,pattern,perl=perl,...)))[1]
406  }
407  if (!is.null(pattern2)) {
408    st <- (unlist(strsplit(st,pattern2,perl=perl,...)))
409    if (reverse) st <- rev(st)
410  }
411  unlist(strsplit(st,sep,perl=perl,...))
412}
413
414Decomp.specials <- function(x,pattern="[()]") {
415  st <- gsub(" ","",x)
416  st <- gsub("\n","",st)
417  mysplit <- rev(unlist(strsplit(st,pattern)))
418  type <- mysplit[2]
419  vars <- mysplit[1]
420  res <- unlist(strsplit(vars,","))
421  if (type=="s" | type=="seq") {
422    return(paste0(res[1],seq(char2num(res[2]))))
423  }
424  unlist(strsplit(vars,","))
425
426}
427
428printline <- function(n=70) {
429    cat(rep("_", n), "\n", sep="");
430
431}
432