# solved by generalized estimation equation #library(pracma) library(neldermead) #'iter' qenv_gee <- function(X, Y, u,tau=0.5, maxiter = 1e2,ftol=1e-2,init=NULL,asy=0,verbose=0,method=1) { t1 = proc.time() X <- as.matrix(X) Y <- as.matrix(Y) n <- nrow(Y) r <- ncol(Y) p <- ncol(X) mX = colMeans(X) mY = colMeans(Y) Yc <- as.matrix(scale(Y, center = T, scale = FALSE)) Xc <- as.matrix(scale(X, center = T, scale = FALSE)) sigX <- cov(Xc) * (n-1)/n W = cbind(rep(1,n),X) ModelOutput=list() if(u == 0){ ModelOutput$alpha = quantile(Y,tau); ModelOutput$beta = matrix(0,p, r) ModelOutput$Gamma = NULL ModelOutput$Gamma0 = diag(p); ModelOutput$eta = NULL; ModelOutput$SigX = sigX ModelOutput$Omega = NULL; ModelOutput$Omega0 = sigX; #ModelOutput$logobj = - n * (r + p) / 2 * (1 + log(2 * pi)) - n / 2 * (logDetsigX + logDetsigY); ModelOutput$paramNum = r + p + p * (p + 1) / 2 Z = Y-ModelOutput$alpha logobj = log(sum(Z*(tau-(Z<0)))/n) BIC = logobj+log(n)*ModelOutput$paramNum/(2*n) ModelOutput$logobj=logobj ModelOutput$BIC=BIC ModelOutput$n = n; if(asy==1){ ModelOutput$covMatrix = NULL; ModelOutput$asySE = NULL; ModelOutput$ratio = matrix(1,p, r); } } else if (u == p){ ModelOutput=list() m1= rq(Y~X,tau = tau) print(alpha <- m1$coefficients[1]) beta = matrix(m1$coefficients[-1],p,1) muX = mX SigmaX = sigX theta0 = c(alpha,muX,beta,SigmaX) # #========= f<-function(theta){ alpha=theta[1] numpar=1 muX=theta[1:p+numpar] numpar=numpar+p beta = matrix(theta[1:p+numpar],p,1) numpar = numpar+p SigmaX = inv.vech(theta[1:(p*(p+1)/2)+numpar]) ones=matrix(1,n,1) SX <- t(X-ones%*%muX)%*%(X-ones%*%muX)/n h1 = t(W)%*%((Y<(X%*%beta+alpha))-tau)/n h2 = muX-mX h3 = vech( SigmaX-SX) c(h1,h2,h3) } #====== #f<-function(theta) sum(h(theta)^2) #m1=fminsearch(f=f,x0=theta0) if(method==1) theta=fsolve(f =f,x0=theta0,tol=ftol)$x #else if(method==2) theta=newtonsys(Ffun =f,x0=theta0,tol=ftol)$x #f(m1$xval) #theta=m0$xval print((alpha<-theta[1])) numpar=1 muX=theta[1:p+numpar] numpar=numpar+p beta = matrix(theta[1:p+numpar],p,1) numpar = numpar+p SigmaX = inv.vech(theta[1:(p*(p+1)/2)+numpar]) numpar=r + p * r + p + p * (p + 1) / 2; Z = Y-(X%*%beta+alpha) logobj = log(sum(Z*(tau-(Z<0)))/n) BIC = logobj+log(n)*numpar/(2*n) ModelOutput$alpha = alpha ModelOutput$beta = beta #ModelOutput$Gamma = Gamma ModelOutput$paramNum = r + p * r +p + p * (p + 1) / 2; #ModelOutput$Gamma0 = Gamma0 #ModelOutput$eta = eta #ModelOutput$SigX = SigX #ModelOutput$iter = iter #ModelOutput$paramNum = r + u * r +p + p * (p + 1) / 2; #ModelOutput$Gamma0 = Gamma0 #ModelOutput$eta = eta #ModelOutput$SigX = SigX #ModelOutput$iter = iter #ModelOutput$Omega = Omega #ModelOutput$Omega0 = Omega0 #ModelOutput$sigYcX = sigYcX ModelOutput$logobj = logobj ModelOutput$BIC = BIC ModelOutput$paramNum = r + p* r +p + p * (p + 1) / 2 ModelOutput$n = n } else{ #initialization if(is.null(init)) { #m1 = rq(Y~X,tau) #betaQ=m1$coefficients[-1] #Gamma_init=get_Init_quantile(sigX, betaQ%*%t(betaQ),u) Gamma_init = unclass(pls::plsr(Y~X,ncomp=u)$loadings) } else Gamma_init = init cat('Angle of init:',subspace(GammaT,Gamma_init),'\n') W = cbind(rep(1,n),X) G_init = Gamma_init %*% solve(Gamma_init[1:u,]) A_init = G_init[-(1:u),] tmp=qr.Q(qr(G_init), complete = TRUE) Gamma_init<- tmp[, 1:u,drop=FALSE] Gamma0_init<- tmp[, (u+1):p,drop=FALSE] m1 = rq(Y~X%*%G_init,tau) alpha <- m1$coefficients[1] eta_ast = matrix(m1$coefficients[-1],u,1) muX = mX SX = sigX Omega = t(Gamma_init) %*% SX %*% Gamma_init Omega0 = t(Gamma0_init) %*% SX %*% Gamma0_init theta0 = c(alpha,muX,A_init,eta_ast,vech(Omega),vech(Omega0)) # #========= f<-function(theta){ alpha=theta[1] numpar=1 muX=theta[1:p+numpar] numpar=numpar+p A = matrix(theta[1:((p-u)*u)+numpar],p-u,u) numpar = numpar + (p-u)*u eta_ast = matrix(theta[1:u+numpar],u,1) numpar = numpar+u Omega = inv.vech(theta[1:(u*(u+1)/2)+numpar]) numpar = numpar+u*(u+1)/2 Omega0 = inv.vech(theta[1:((p-u)*(p-u+1)/2)+numpar]) G = rbind(diag(u),A) tmp=qr.Q(qr(G), complete = TRUE) Gamma <- tmp[, 1:u,drop=FALSE] Gamma0 <- tmp[, (u+1):p,drop=FALSE] ones=matrix(1,n,1) SX <- t(X-ones%*%muX)%*%(X-ones%*%muX)/n h1 = t(W)%*%((Y<(X%*%G%*%eta_ast+alpha))-tau)/n h2 = muX-mX h3 = vech( Gamma %*%Omega%*%t(Gamma)+Gamma0%*% Omega0 %*%t(Gamma0)-SX) c(h1,h2,h3) } #====== #f<-function(theta) sum(h(theta)^2) #m1=fminsearch(f=f,x0=theta0) if(method==1) theta=fsolve(f =f,x0=theta0,tol=ftol)$x alpha<-theta[1] numpar=1 muX=theta[1:p+numpar] numpar=numpar+p A = matrix(theta[1:((p-u)*u)+numpar],p-u,u) numpar = numpar + (p-u)*u eta_ast = matrix(theta[1:u+numpar],u,1) numpar = numpar+u Omega = inv.vech(theta[1:(u*(u+1)/2)+numpar]) numpar = numpar+u*(u+1)/2 Omega0 = inv.vech(theta[1:((p-u)*(p-u+1)/2)+numpar]) numpar = numpar+(p-u)*(p-u+1)/2 G = rbind(diag(u),A) #print(G%*%eta_ast-betaT) tmp=qr.Q(qr(G), complete = TRUE) Gamma <- tmp[, 1:u,drop=FALSE] Gamma0 <- tmp[, (u+1):p,drop=FALSE] ones=matrix(1,n,1) SX <- t(X-ones%*%muX)%*%(X-ones%*%muX)/n cat('Angle of gee estimator:',subspace(GammaT,Gamma),'\n') #G = G_init numpar=r + u * r + p + p * (p + 1) / 2; Z = Y-(X%*%G%*%eta_ast+alpha) logobj = log(sum(Z*(tau-(Z<0)))/n) BIC = logobj+log(n)*numpar/(2*n) ModelOutput$alpha = alpha ModelOutput$beta = G%*%eta_ast ModelOutput$Gamma = Gamma ModelOutput$paramNum = r + u * r +p + p * (p + 1) / 2; #ModelOutput$Gamma0 = Gamma0 #ModelOutput$eta = eta #ModelOutput$SigX = SigX #ModelOutput$iter = iter #ModelOutput$Omega = Omega #ModelOutput$Omega0 = Omega0 #ModelOutput$sigYcX = sigYcX ModelOutput$logobj = logobj ModelOutput$BIC = BIC #ModelOutput$paramNum = r + u * r +p + p * (p + 1) / 2 ModelOutput$n = n if(asy==1){ #---compute asymptotic variance and get the ratios--- asyFm = kronecker(sigYcX, invsigX); asyFm = matrix(sqrt(diag(asyFm)), p, r); temp = kronecker(eta %*% solve(sigYcX) %*% t(eta), Omega0)+ kronecker(Omega, solve(Omega0)) + kronecker(solve(Omega), Omega0) - 2 * kronecker(diag(u), diag(p - u)) covMatrix = kronecker(sigYcX, Gamma %*% solve(Omega) %*% t(Gamma)) + kronecker(t(eta), Gamma0) %*% solve(temp) %*% kronecker(eta, t(Gamma0)); asySE = matrix(sqrt(diag(covMatrix)), p, r) ModelOutput$covMatrix = covMatrix ModelOutput$asySE = asySE ModelOutput$ratio = asyFm / asySE } } class(ModelOutput) <- "qenv_gee" ModelOutput$fit_time=(proc.time()-t1)[3] return(ModelOutput) }