Ds <- list(structure(c(3, 2, 2, 4), .Dim = c(2L, 2L)), structure(c(1, 0, 0, 2), .Dim = c(2L, 2L)), structure(c(3, 0, 0, 3), .Dim = c(2L, 2L))) lmm <- function (Y0, Zt, Zs0, covType, PRINT=F, IV=T, Ds=NULL, s2=NULL, PLOT=F, MAXITR=1000){ L = nrow(Y0) sh = rep(1,L) alpha = apply(Y0*sh,2,mean) Y = t(t(Y0)-alpha)*sqrt(sh) print(dim(Y)) Vh = t(Y)%*%Y/(L-1) tr = function(x){sum(diag(x))} asum = function(x){sum(abs(x))} Zs0 = lapply(Zs0,function(z){matrix(z[,abs(apply(z,2,sum))>0],nrow(z))}) #Zt = matrix(Zt[,apply(matrix(Zt,46),2,sum)>0],nrow(matrix(Zt,46))) K = ncol(Zt) N = ncol(Y) Zall=NULL; for(i in c(1:length(Zs0))){Zall=cbind(Zall, Zs0[[i]])} Q=qr.Q(qr(Zall)) if(is.null(s2)){ ######### INIT s2 s2 = NULL m=apply(matrix(Y,L),1,mean) for(i in 1:K){ s2 = c(s2, var(c(Y[,Zt[,i]==1]-m))) } #s2=s20 if(PRINT){print(s2)} T1=1/sqrt(c(Zt%*%s2)) } ######## INIT Zs Zs = Zs0 for(i in 1:length(Zs)){ #Zs[[i]] = Zs0[[i]]*T1 } if(is.null(Ds)){ ########## INIT Ds Ds=NULL for(i in 1:length(Zs)){ Yi = NULL for(k in 1:ncol(Zs[[i]])){ Yi = cbind(Yi, apply(matrix(Y[,Zs[[i]][,k]>0],L),1,mean)) } if(covType[i]%in%c(1,4)){ Ds = c(Ds, list(var(Yi)/length(Zs))) }else if(covType[i]==2){ Ds = c(Ds, list(diag(diag(var(Yi)/length(Zs))))) }else if(covType[i]==2.5||covType[i]==3.5||covType[i]==3){ Ds = c(Ds, list(diag(ncol(Yi))*mean(diag(var(Yi)/length(Zs))))) }else if(covType[i]==1.5){ Ds = c(Ds, list(Dg*0.4)) } for(k in 1:ncol(Zs[[i]])){ #Y[,Zs[[i]][,k]==1] = Y[,Zs[[i]][,k]==1] - apply(matrix(Y[,Zs[[i]][,k]==1],L),1,mean) } } #diag(Ds[[2]])=abs(diag(Ds[[2]])-c(Ds[[1]])) for(i in 2:length(Ds)){Ds[[i]]=Ds[[i]]/100} } #Ds[[10]]=Ds[[9]]=matrix(0.001,1,1) if(PRINT){print(Ds)} NU = 1.526 ######################### if(PRINT){cat("iteration......\n")} lkhd=NULL;Us=NULL;Vt=NULL;V=NULL;T2=NULL for (itr in 1:MAXITR) { cat(paste("Iteration: ",itr," Likelihood: ",lkhd[length(lkhd)],"\n",sep="")) ngrad=0 if(IV){ T2 = 1/c(Zt%*%s2) for(i in 1:length(Zs)){Zs[[i]] = Zs0[[i]]*sqrt(T2)} Vt = diag(N); for(l in 1:length(Zs)){Vt = Vt + Zs[[l]] %*% Ds[[l]] %*% t(Zs[[l]])} V = t(Vt/sqrt(T2))/sqrt(T2) Vinv=Solve(V) alpha = apply(Y0*sh,2,sum)/sum(sh) RR=try(chol(V)) if(is.character(RR)){RR=chol((V+diag(diag(V)))/2)} Q1=forwardsolve(t(RR),Q) alpha=forwardsolve(t(RR),alpha) #print(dim(Q)) #print(as.numeric(coef(lm(alpha~Q1)))) beta=c(Q%*%as.numeric(coef(lm(alpha~Q1-1)))) #beta=alpha sh = apply(t(t(Y0)-beta),1,function(yh){sum(t(Vinv*yh)*yh)}) lsh = digamma((ncol(Y)+NU)/2)-log((sh+NU)/2) sh = (ncol(Y)+NU)/(sh+NU) msh=mean(sh) mlsh=mean(lsh) #sh = sh/mean(sh) #alpha = apply(Y0*sh,2,sum)/sum(sh) #RR=try(chol(V)) #if(is.character(RR)){RR=chol((V+diag(diag(V)))/2)} #Q1=forwardsolve(t(RR),Q) #alpha=forwardsolve(t(RR),alpha) #beta=c(t(RR)%*%predict(lm(alpha~Q1))) #beta=beta*0 Y = t(t(Y0)-beta)*sqrt(sh) Vh = t(Y)%*%Y/(L-1) if(itr>10){ for(itrNU in 1:1000){ g = (mlsh-msh+log(NU/2)+1-digamma(NU/2))/2 h = (1/NU-trigamma(NU/2)/2)/2 if(g>0){ NU = exp( log(NU) + abs(g/(h*NU+g)) ) }else{ NU = exp( log(NU) - abs(g/(h*NU+g)) ) } #if(PRINT)print(NU) } } if(NU<0.1 || is.na(NU)){NU=2.21872} } #print(Ds[[1]]) if(Ds[[2]][1]>1.4){diag(Ds[[2]])=0.6} #if(Ds[[1]][1]100){ for (itrU in 1:500) { Vt = diag(N); for(l in 1:length(Zs)){Vt = Vt + Zs[[l]] %*% Ds[[l]] %*% t(Zs[[l]])} V = t(Vt/sqrt(T2))/sqrt(T2) A = t(Zs[[k]]) %*% solve(Vt, Zs[[k]]) #B = t(Zs[[k]]/sqrt(T2)) %*% solve(V, Vh) %*% solve(sqrt(T2)*V, Zs[[k]]) B = t(Zs0[[k]]) %*% solve(V, Vh) %*% solve(V, Zs0[[k]]) U = Eigen(Ds[[k]]) g = c(t((B - A) %*% t(U))) h = A %x% (U %*% A %*% t(U)) + (A %*% t(U)) %x% (U %*% A) %*% getK(ncol(Zs[[k]])) dU = c(Solve(h) %*% c(g)) lkhd0 = -log(det(V))-sum(diag(solve(V,Vh))) for(steps in 0:(-10)){ U1 = U + dU*(2^steps) Ds[[k]] = t(U1) %*% U1 Vt = diag(N); for(l in 1:length(Zs)){Vt = Vt + Zs[[l]] %*% Ds[[l]] %*% t(Zs[[l]])} V = t(Vt/sqrt(T2))/sqrt(T2) #if(is.na(det(V))){print(Ds[[k]])} if(lkhd0<(-log(det(V))-sum(diag(solve(V,Vh))))){break} } U = U1 Ds[[k]] = t(U) %*% U if(asum(g)<1e-8){ses=c(ses,list(diag(Solve(h))/(L-1))); break} lkhd1=c(lkhd1,(-log(det(V))-sum(diag(solve(V,Vh))))) if(itrU>1){if(abs(lkhd1[length(lkhd1)]-lkhd1[length(lkhd1)-1])<1e-10){ses=c(ses,list(diag(Solve(h))/(L-1))); ngrad=ngrad+1; break}} } if(PRINT){cat(1);print(asum(g))} #}else{U=Eigen(Ds[[k]]);Ds[[k]]=matrix(c(20,13,13,14),2)} }else if(covType[k]==1.5){ for (itrU in 1:500) { Vt = diag(N); for(l in 1:length(Zs)){Vt = Vt + Zs[[l]] %*% Ds[[l]] %*% t(Zs[[l]])} V = t(Vt/sqrt(T2))/sqrt(T2) A = t(Zs[[k]]) %*% solve(Vt, Zs[[k]]) B = t(Zs0[[k]]) %*% solve(V, Vh) %*% solve(V, Zs0[[k]]) U = sqrt(mean(Ds[[k]]/Dg)) ######################## # tr(AB) = vec(A^T)^T vec(B) # U = sigma * R; R = chol(D) # g = vec( sigma * ( (B-A) R^T )^T )^T vec(R) # = sum(sigma * diag((B-A)D)) # h = tr(V^-1 dV V^-1 dV) # = tr( V^-1 Z D Z^T V^-1 Z D Z^T ) * sigma^2 # = tr( (Z^T V^-1 Z D)^2 ) * sigma^2 ######################## g = sum(diag((B - A) %*% Dg))*U h = sum(diag( (A%*%Dg)%*%(A%*%Dg) ))*U^2 if(g*U>0){dU = abs(g*U/(h*U^2+g*U))}else{dU = -abs(g*U/(h*U^2+g*U))} lkhd0 = -log(det(V))-sum(diag(solve(V,Vh))) for(steps in 0:(-10)){ U1 = exp(log(U)+dU*(2^steps)) Ds[[k]] = Dg * U1^2 Vt = diag(N); for(l in 1:length(Zs)){Vt = Vt + Zs[[l]] %*% Ds[[l]] %*% t(Zs[[l]])} V = t(Vt/sqrt(T2))/sqrt(T2) #if(is.na(det(V))){print(Ds[[k]])} if(lkhd0<(-log(det(V))-sum(diag(solve(V,Vh))))){break} } U = U1 Ds[[k]] = Dg * U if(asum(g)<1e-8){ses=c(ses,list(1/sum(h*(L-1)))); break} lkhd1=c(lkhd1,(-log(det(V))-sum(diag(solve(V,Vh))))) if(itrU>1){if(abs(lkhd1[length(lkhd1)]-lkhd1[length(lkhd1)-1])<1e-10){ses=c(ses,list(1/sum(h*(L-1)))); ngrad=ngrad+1; break}} } if(PRINT){cat(1);print(asum(g))} }else if(covType[k]==2){# diagonal different var #if(PRINT){cat("Batch ");print(Ds[[k]])} for (itrU in 1:500) { M = ncol(Zs[[k]]) Vt = diag(N); for(l in 1:length(Zs)){Vt = Vt + Zs[[l]] %*% Ds[[l]] %*% t(Zs[[l]])} V = t(Vt/sqrt(T2))/sqrt(T2) A = t(Zs0[[k]]) %*% solve(V, Zs0[[k]]) B = t(Zs0[[k]]) %*% solve(V, Vh) %*% solve(V, Zs0[[k]]) U = sqrt(Ds[[k]]) ######################## # U = sigma * I # g_i = vec( ( (B-A) diag(sigma_i; i=1,...)^T )^T )^T vec(E_i) ; E_i : (e_ii=1 and other=0) # = diag(B-A)*(sigma_i) ######################## g = diag(t((B - A) %*% t(U))) h = (A %x% (U %*% A %*% t(U)) + (A %*% t(U)) %x% (U %*% A) %*% getK(M))[c(diag(M))==1,c(diag(M))==1] dU = c(Solve(h) %*% c(g)) lkhd0 = -log(det(V))-sum(diag(solve(V,Vh))) for(steps in 0:(-10)){ U1 = U + diag(dU*(2^steps)) Ds[[k]] = U1^2 Vt = diag(N); for(l in 1:length(Zs)){Vt = Vt + Zs[[l]] %*% Ds[[l]] %*% t(Zs[[l]])} V = t(Vt/sqrt(T2))/sqrt(T2) #if(is.na(log(det(V)))){print(Ds[[k]]);print(k)} if(lkhd0<(-log(det(V))-sum(diag(solve(V,Vh))))){break} } U = U1 Ds[[k]] = U^2 if(asum(g)<1e-8){ses=c(ses,list(diag(Solve(h))/(L-1)));break} lkhd1=c(lkhd1,(-log(det(V))-sum(diag(solve(V,Vh))))) if(itrU>1){if(abs(lkhd1[length(lkhd1)]-lkhd1[length(lkhd1)-1])<1e-10){ses=c(ses,list(diag(Solve(h))/(L-1))); ngrad=ngrad+1; break}} } if(PRINT){cat(2);print(asum(g))} #if(PRINT){print(Ds[[k]])} }else if(covType[k]==3){# diagonal same var #if(PRINT){cat("Batch ");print(Ds[[k]])} for (itrU in 1:500) { M = ncol(Zs[[k]]) Vt = diag(N); for(l in 1:length(Zs)){Vt = Vt + Zs[[l]] %*% Ds[[l]] %*% t(Zs[[l]])} V = t(Vt/sqrt(T2))/sqrt(T2) A = t(Zs0[[k]]) %*% solve(V, Zs0[[k]]) A = (A+t(A))/2 B = t(Zs0[[k]]) %*% solve(V, Vh) %*% solve(V, Zs0[[k]]) U = sqrt(Ds[[k]]) #g = diag(t((B - A) %*% t(U))) ######################## # U = sigma * I # g = vec( ( sigma * (B-A) I^T )^T )^T vec(I) # = sum(diag(B-A)*sigma) ######################## g = diag(B-A)*U[1] #h = (A %x% (U %*% A %*% t(U)) + (A %*% t(U)) %x% (U %*% A) %*% getK(M))[c(diag(M))==1,c(diag(M))==1] h = diag(2*A%*%A*(U[1]^2)) if(sum(g)*U[1]>0){dU = abs(sum(g)*U[1]/(sum(h)*U[1]^2+sum(g)*U[1]))}else{dU = -abs(sum(g)*U[1]/(sum(h)*U[1]^2+sum(g)*U[1]))} lkhd0 = -log(det(V))-sum(diag(solve(V,Vh))) for(steps in 0:(-2)){ #U1 = U + diag(nrow(U))*dU*(2^steps) U1 = diag(nrow(U))*exp(log(U[1])+dU*(2^steps)) Ds[[k]] = U1^2 Vt = diag(N); for(l in 1:length(Zs)){Vt = Vt + Zs[[l]] %*% Ds[[l]] %*% t(Zs[[l]])} V = t(Vt/sqrt(T2))/sqrt(T2) #if(is.na(log(det(V)))){print(Ds[[k]]);print(k)} if(lkhd0<(-log(det(V))-sum(diag(solve(V,Vh))))){break} } U = U1 Ds[[k]] = U^2 if(sum(g)*U[1]<1e-8){ses=c(ses,list(1/sum(h*(L-1))));break} lkhd1=c(lkhd1,(-log(det(V))-sum(diag(solve(V,Vh))))) #if(itrU>1){if(abs(lkhd1[length(lkhd1)]-lkhd1[length(lkhd1)-1])<1e-18){ses=c(ses,list(1/sum(h*(L-1)))); ngrad=ngrad+1; break}} } if(PRINT){cat(2.5);print(sum(g))} #if(PRINT){print(Ds[[k]])} }else if(covType[k]==3.5){# diagonal different var partitioned #if(PRINT){cat("Batch ");print(Ds[[k]])} if(k==3){ Zz=diag(2)[c(1,1,1,2,2),] }else if(k==4){ Zz=diag(2)[c(1,1,2),] }else if(k==5){ Zz=diag(2)[c(1,1,1,2,2,2,2,2),] }else if(k==6){ Zz=Zz.brep2#[ brep[30,]==0,] } for (itrU in 1:500) { M = ncol(Zs[[k]]) Vt = diag(N); for(l in 1:length(Zs)){Vt = Vt + Zs[[l]] %*% Ds[[l]] %*% t(Zs[[l]])} V = t(Vt/sqrt(T2))/sqrt(T2) A = t(Zs0[[k]]) %*% solve(V, Zs0[[k]]) B = t(Zs0[[k]]) %*% solve(V, Vh) %*% solve(V, Zs0[[k]]) U = sqrt(Ds[[k]]) g = diag(t((B - A) %*% t(U))) h = (A %x% (U %*% A %*% t(U)) + (A %*% t(U)) %x% (U %*% A) %*% getK(M))[c(diag(M))==1,c(diag(M))==1] dU = Solve(t(Zz)%*%h%*%Zz)%*%t(Zz)%*%g lkhd0 = -log(det(V))-sum(diag(solve(V,Vh))) for(steps in 0:(-10)){ U1 = U + diag(c(Zz%*%dU))*(2^steps) Ds[[k]] = U1^2 Vt = diag(N); for(l in 1:length(Zs)){Vt = Vt + Zs[[l]] %*% Ds[[l]] %*% t(Zs[[l]])} V = t(Vt/sqrt(T2))/sqrt(T2) #if(is.na(log(det(V)))){print(Ds[[k]]);print(k)} if(lkhd0<(-log(det(V))-sum(diag(solve(V,Vh))))){break} } U = U1 Ds[[k]] = U^2 if(asum(t(Zz)%*%g)<1e-8){ses=c(ses,list(diag(Solve(t(Zz)%*%h%*%Zz))/(L-1))); break} lkhd1=c(lkhd1,(-log(det(V))-sum(diag(solve(V,Vh))))) if(itrU>1){if(abs(lkhd1[length(lkhd1)]-lkhd1[length(lkhd1)-1])<1e-10){ses=c(ses,list(diag(Solve(t(Zz)%*%h%*%Zz))/(L-1))); ngrad=ngrad+1; break}} } if(PRINT){cat(3.5);print(asum(g))} #if(PRINT){print(Ds[[k]])} }else if(covType[k]==3.6){# diagonal different var patitioned with correlation R Zz=diag(2)[c(1,1,1,2,2,2,2,2),] for (itrU in 1:500) {# iteration for S M = ncol(Ds[[k]]) S = sqrt(diag(Ds[[k]])) U = S Vt = diag(N); for(l in 1:length(Zs)){Vt = Vt + Zs[[l]] %*% Ds[[l]] %*% t(Zs[[l]])} V = t(Vt/sqrt(T2))/sqrt(T2) A = t(Zs0[[k]]) %*% solve(V, Zs0[[k]]) B = t(Zs0[[k]]) %*% solve(V, Vh) %*% solve(V, Zs0[[k]]) SDA = (Ds[[k]]/S)%*%A g = diag((B - A) %*% t(Ds[[k]]/S)) h = (A%x%(SDA%*%t(Ds[[k]]/S)) + (SDA%x%t(SDA))%*%getK(M) )[c(diag(M))==1,c(diag(M))==1] dU = Solve(t(Zz)%*%h%*%Zz)%*%t(Zz)%*%g lkhd0 = -log(det(V))-sum(diag(solve(V,Vh))) for(steps in 0:(-10)){ U1 = abs(U + c(Zz%*%dU)*(2^steps)) Ds[[k]] = t(R*U1)*U1 Vt = diag(N); for(l in 1:length(Zs)){Vt = Vt + Zs[[l]] %*% Ds[[l]] %*% t(Zs[[l]])} V = t(Vt/sqrt(T2))/sqrt(T2) if(lkhd0<(-log(det(V))-sum(diag(solve(V,Vh))))){break} } U = U1 Ds[[k]] = t(R*U)*U if(asum(g)<1e-8){ses=c(ses,list(diag(Solve(h))/(L-1))); break} lkhd1=c(lkhd1,(-log(det(V))-sum(diag(solve(V,Vh))))) if(itrU>1){if(abs(lkhd1[length(lkhd1)]-lkhd1[length(lkhd1)-1])<1e-10){ses=c(ses,list(diag(Solve(h))/(L-1))); ngrad=ngrad+1; break}} } #if(PRINT){print(Cor(Ds[[k]]))} if(PRINT){cat(3);print(asum(g))} } Us = c(Us, list(U)) #plot(lkhd1) #print(c(itrU,sum(abs(g)))) } # iteration for T for(i in 1:length(Zs)){Zs[[i]] = Zs0[[i]]*sqrt(T2)} Vt = diag(N); for(l in 1:length(Zs)){Vt = Vt + Zs[[l]] %*% Ds[[l]] %*% t(Zs[[l]])} V = t(Vt/sqrt(T2))/sqrt(T2) for(itrT in 1:1000){ T2 = 1/c(Zt%*%s2) g = apply(Zt*sqrt(1/T2),2,sum) - apply(Zt*diag(solve(Vt, Vh*sqrt(T2))), 2, sum) tmp = - Vh*solve(Vt) - diag(1/T2) h = array(0, c(K,K)) for(k in 1:K){for(l in 1:K){h[k,l] = sum(tmp[Zt[,k]==1,Zt[,l]==1])}} s2 = c(1/(1/sqrt(s2) - Solve(h)%*%g)^2) if(sum(abs(g))<1e-8){ses=c(ses,list(-diag(Solve(h))/(L-1))); break} } #print("tau") #print(g) #print(c(itrT,sum(abs(g)))) #print("Ts");print(sum(abs(g))) T2 = 1/c(Zt%*%s2) for(i in 1:length(Zs)){Zs[[i]] = Zs0[[i]]*sqrt(T2)} Vt = diag(N); for(l in 1:length(Zs)){Vt = Vt + Zs[[l]] %*% Ds[[l]] %*% t(Zs[[l]])} V = t(Vt/sqrt(T2))/sqrt(T2) lkhd=c(lkhd,(L-1)/2*(-log(det(V))-sum(diag(solve(V,Vh))))) #print(lkhd[length(lkhd)]) if(PLOT)plot(lkhd) #barplot(c(getDiag(Ds[-1]),s2)) if(itr>2){if(abs(lkhd[length(lkhd)]-lkhd[length(lkhd)-1])<1e-8){break}} hoge = list(s2, Ds, list(V,Vh), Vt, ses, sh) #if(itr%%100==0){dump("hoge","hoge.R")} #print(sum(abs(g))) } Vt = list(diag(N)/T2); for(l in 1:length(Zs)){Vt = c(Vt , list(Zs0[[l]] %*% Ds[[l]] %*% t(Zs0[[l]])))} #print(lkhd[length(lkhd)]) colnames(V)=colnames(Vh) rownames(V)=rownames(Vh) list(s2, Ds, list(V,Vh), Vt, ses, sh, NU) } Sigma <- structure(c(1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 2), .Dim = c(8L, 8L)) X <- structure(list(x1 = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L ), .Label = c("1", "2"), class = "factor"), x2 = structure(c(1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L), .Label = c("1", "2"), class = "factor"), x3 = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("1", "2"), class = "factor")), .Names = c("x1", "x2", "x3"), row.names = c(NA, 8L), class = "data.frame") Solve <- function(x,th=1e-9){ if(nrow(x)==1){return(1/x)} x=eigen((x+t(x))/2,T); x[[1]][abs(x[[1]])th]=1/x[[1]][abs(x[[1]])>th]; x[[2]]%*%diag(x[[1]])%*%t(x[[2]]) } Eigen <- function(x){ if(nrow(x)==1){return(sqrt(x))} x=(x+t(x))/2 x=eigen(x,T); x[[1]][x[[1]]<0]=0; diag(sqrt(x[[1]]))%*%t(x[[2]])} getDiag <- function(x){x=lapply(x,function(xx){c(unique(diag(xx)),NA)});unlist(x)} getK <- function(n){diag(n^2)[c(t(matrix(1:(n^2),n))),]}