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