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