1## Make random effects for mixed-effect models 2mkran <- function(formula,data) 3{ 4 ## decipher formula 5 form.wk <- terms.formula(formula)[[2]] 6 terms <- strsplit(deparse(form.wk),' \\+ ')[[1]] 7 if (length(terms)>1) { 8 form <- as.formula(paste("~",terms[1])) 9 zzz <- mkran(form,data) 10 for (i in 2:length(terms)) { 11 form <- as.formula(paste("~",terms[i])) 12 zzz <- mkran1(zzz,mkran(form,data)) 13 } 14 return(zzz) 15 } 16 if (!("|"%in%strsplit(deparse(form.wk),'')[[1]])) 17 stop("gss error in mkran: missing | in grouping formula") 18 term.wk <- strsplit(deparse(form.wk),' \\| ')[[1]] 19 with(data,{ 20 ## make matrix Z 21 z2.wk <- eval(parse(text=term.wk[2])) 22 if (!is.factor(z2.wk)) 23 stop(paste("gss error in mkran: ", term.wk[2], " should be a factor")) 24 z <- NULL 25 lvl.z2 <- levels(z2.wk) 26 for (i in lvl.z2) z <- cbind(z,as.numeric(z2.wk==i)) 27 ## make sigma function 28 if (term.wk[1]=="1") { 29 init <- 0 30 env <- length(levels(z2.wk)) 31 fun <- function(zeta,env) diag(10^(-zeta),env) 32 sigma <- list(fun=fun,env=env) 33 } 34 else { 35 z1.wk <- eval(parse(text=term.wk[1])) 36 if (!is.factor(z1.wk)) 37 stop(paste("gss error in mkran: ", term.wk[1], " should be a factor")) 38 ind <- lvl.wk <- NULL 39 nz <- length(lvl.z2) 40 nsig <- length(levels(z1.wk)) 41 for (i in levels(z1.wk)) { 42 zz.wk <- z2.wk[z1.wk==i,drop=TRUE] 43 ind <- c(ind,list((1:nz)[lvl.z2%in%levels(zz.wk)])) 44 lvl.wk <- c(lvl.wk,levels(zz.wk)) 45 } 46 if (max(table(lvl.wk)>1)) 47 stop("gss error in mkran: ", term.wk[2], " should be nested under ", term.wk[1]) 48 init <- rep(0, length(levels(z1.wk))) 49 env <- list(size=nz,nsig=nsig,ind=ind) 50 fun <- function(zeta,env) { 51 wk <- rep(0,env$size) 52 for (i in 1:env$nsig) wk[env$ind[[i]]] <- 10^(-zeta[i]) 53 diag(wk) 54 } 55 sigma <- list(fun=fun,env=env) 56 } 57 list(z=z,sigma=sigma,init=init) 58 }) 59} 60 61## Combine random effects for mixed-effect models 62mkran1 <- function(ran1,ran2) 63{ 64 z <- cbind(ran1$z,ran2$z) 65 env <- list(sz1=dim(ran1$z)[2],sig1=ran1$sigma,nz1=length(ran1$init), 66 sz2=dim(ran2$z)[2],sig2=ran2$sigma,nz2=length(ran2$init)) 67 fun <- function(zeta,env) { 68 idx1 <- 1:env$sz1 69 idx2 <- env$sz1+(1:env$sz2) 70 sig <- matrix(0,env$sz1+env$sz2,env$sz1+env$sz2) 71 sig[idx1,idx1] <- env$sig1$fun(zeta[1:env$nz1],env$sig1$env) 72 sig[idx2,idx2] <- env$sig2$fun(zeta[env$nz1+(1:env$nz2)],env$sig2$env) 73 sig 74 } 75 sigma <- list(fun=fun,env=env) 76 list(z=z,sigma=sigma,init=c(ran1$init,ran2$init)) 77} 78