R/CLcov.R

In m-azad/ComLikMulCom: Multiple comparisons using composite likelihood

CLcov <- function(data, yno , idno, N, n, p , y.type = c("normal", "probit", "quadratic exponential")[1] ,  type = c("many_one", "all_pair")[1]){
  
  require(mvtnorm)
  
  data <- data.frame(data)
  
  yy <- as.matrix(data[,yno])
  xx <- as.matrix(data[,-c(yno,idno)])
  
  

if("normal" %in% y.type){
  ###    Estimating beta and sigma from marginal normal log likelihood
  
  sum1     <-   t(xx)%*%xx
  sum2     <-   t(xx)%*%yy
  beta_hat <-   solve(sum1)%*%sum2  
  
  
  yhat      <-   matrix(0,nrow=m,ncol=n)
  e         <-   matrix(0,nrow=m,ncol=n)
  
  for(i in 1: n){
    index <- which(data[, idno]==i)
    d <- as.matrix(data[index,])
    e[,i] <- d[,yno] - d[,-c(idno,yno)]%*%beta_hat
  }
  
  
  sigma_hat <-   matrix(0, nrow=m, ncol=m)
  
  
  sii     <-   sum(apply(e,1,function(x)sum(x^2)))/(n*m-p)	
  
  E         <- 	rep(0,n)
  for (i in 1:n){
    for( j in 1:(m-1)){
      E[i] 	<- 	E[i] + sum(e[j,i]*e[(j+1):m,i])
    }  
  }
  r 		<- 	sum(E)/((n-p)*choose(m,2)*sii)
  
  Sigma_hat  <- 	array(rep(r,m^2),dim=c(m,m))
  Sigma_hat  <- 	Sigma_hat+ (1-r)*diag(m)
  Sigma_hat  <- 	sii*Sigma_hat 		                 
  
  ###  	computing covariance of coefficients
  
  c1		<-	matrix(0,nrow=p,ncol=p)
  c2		<- 	matrix(0,nrow=p,ncol=p)
  
  for(i in 1:n){
    
    index <- which(data[, idno]==i)
    d <- as.matrix(data[index,])
    
    #d <- as.matrix(data[data$id==i,])
    c1	<-	c1 +(t(d[,-c(idno,yno)])%*%d[,-c(idno,yno)])
    c2	<-	c2 +t(d[,-c(idno,yno)])%*%Sigma_hat%*%d[,-c(idno,yno)]
  }
  
  Cov_beta  <- 	solve(c1)%*%c2 %*%t(solve(c1)) 
  
  
} 



if("probit" %in% y.type){
  
  fun1_b <- function(s,x,y){  
    
    ps   <- as.numeric(pnorm(s,0,1,lower.tail=TRUE))
    pms  <- as.numeric(pnorm(-s,0,1,lower.tail=TRUE))
    ds   <- as.numeric(dnorm(s,0,1))  
    if  (abs(s)<= 5) {z1 <- ((y-ps)*(ds*x))/(ps*pms)
    } else if (s>5 & y==1) {z1 <- ds*x/ps
    } else if (s>5 & y==0) {z1 <- -s*x
    } else if (s< -5 & y==1) {z1 <- -s*x
    } else if (s< -5 & y==0) {z1 <- -ds*x/pms}
    return(z1)
  }
  
  
  fun2_b <- function(s,x,y){  
    
    ps   <- as.numeric(pnorm(s,0,1,lower.tail=TRUE))
    pms  <- as.numeric(pnorm(-s,0,1,lower.tail=TRUE))
    ds   <- as.numeric(dnorm(s,0,1))
    if  (y==1 & s>-5 ) {z2 <-as.numeric((-ds^2/ps^2) - (ds*s/ps))*x%*%t(x)
    }else if (y==1 & s< -5) {z2 <- 0
    }else if (y==0 & s< 5 ) {z2 <-as.numeric((-ds^2/pms^2) + (ds*s/pms))*x%*%t(x)
    }else if (y==0 & s> 5) {z2 <- 0}
    return(z2)
  }
  ###Initial value
  
  fit       <-    glm(yy~ xx ,family=binomial(link = "probit"), data , 
                      control = list(epsilon = 1e-6, maxit = 30, trace = FALSE))
  
  beta_hat  <-    as.vector(fit$coefficients[-1])
  
  
  ###Newton Raphson for beta
  
  sig     <-   1
  
  for (k in 1:1000){ 
    
    beta_hat_n   <-  beta_hat 
    sc_b         <-  rep(0,p)
    hs_b         <-  matrix(0,p,p)
    
    for(i in 1:nrow(xx)){
      
      X    <-  xx[i,]%*%beta_hat_n/sig
      
      sc_b <- sc_b + fun1_b(X,xx[i,],yy[i,])
      hs_b <- hs_b + fun2_b(X,xx[i,],yy[i,])
    }
    
    beta_hat <- ifelse(sig*abs(solve(hs_b)%*%sc_b) < 10^(-5),beta_hat_n,beta_hat_n - sig*solve(hs_b)%*%sc_b)
    
    if (all(beta_hat== beta_hat_n)){break}
  } 

  v_s  <- matrix(0,p,p)
  H    <- matrix(0,p,p)
  
  
  for ( i in 1:n){
    
    index <- which(data[, idno]==i)
    d <- as.matrix(data[index,])
    
    dmu   <-   as.vector(dnorm(d[,-c(idno,yno)]%*%beta_hat))*d[,-c(idno,yno)]
    mu    <-   pnorm(d[,-c(idno,yno)]%*%beta_hat)
    V     <-   diag(as.vector(pnorm(d[,-c(idno,yno)]%*%beta_hat))*as.vector(pnorm(-d[,-c(idno,yno)]%*%beta_hat)))
    sij   <-   (d[,yno]-mu)%*%t(d[,yno]-mu)
    
    v_s   <-   v_s + t(dmu)%*%solve(V)%*%sij%*%solve(V)%*%dmu
    H     <-   H + t(dmu)%*%solve(V)%*%dmu
    
  }
  
  Cov_beta  <-  solve(H)%*%v_s%*%solve(H)
  
}

if("quadratic exponential" %in% y.type){
  ## creating the column of dependence factor (w)
  
  pat      <-  array(0,dim=c(n,2))  
  pat[,2]  <-  table(data[,idno])
  ct       <-  1
  breakdown<-  vector("list", n)
  
  
  for (i in 1:(N-1)){
    breakdown[[ct]]  <- c(breakdown[[ct]],data[i,yno])
    if (data[,idno][i+1]>data[,idno][i]){ct<-ct+1}
  }
  for (i in N:N){ 
    breakdown[[ct]]  <- c(breakdown[[ct]],data[i,yno])
  }
  for ( i in 1: n){
    brk <- breakdown[[i]]
    brk <- unlist(brk)
    pat[i,1]         <- sum(brk)
  }
  
  
  ss     <-  pat[,1]
  ff     <-  pat[,2]
  data$z <-  rep(ss,ff)
  data$m <-  rep(ff,ff)
  data$w <-  ifelse(data[,yno]==1,  -(data$m-2*data$z+1),-(data$m-2*data$z-1))
  
  
  r     <-  dim(data)[2]
  xdat  <-  data[,-c((r-1),(r-2))] 
  
  factnam  <-  names(xdat[,-c(idno,yno)])
  
  resp     <-  xdat[,yno]  
  
  fmla     <-  as.formula(paste("resp ~ ", paste(factnam, collapse= "+")))
  fit      <-  glm(fmla ,data=xdat ,family=binomial(link=logit))
  beta <-  fit$coefficients[-1]
  
  beta_hat     <-   as.matrix(fit$coefficients[2:(p+1)])
  what     <-   fit$coefficients[(p+2)]
  
  
  X  <- split(data, data[,idno])
  X  <- lapply(X, function(x) { x[c("id")] <- NULL; x })
  X  <- lapply(X, function(x) { x[c("w")] <- NULL; x })
  
  
  scorevec  <-    array(0,dim=c(n,p))
  
  for ( i in 1: n){
    
    xxs      <-    as.data.frame(X[[i]])
    xxs$indc <-    ifelse(xxs[,yno]==1, 1, -1)
    e        <-    rep(0,nrow(xxs))
    no3      <-    which( colnames(xxs)=="indc" )
    
    for(q in 1: nrow(xxs)){
      med <- xxs[,no3]
      e[q] <- ifelse(med[q]==1,exp(as.matrix(xxs[q,-c(yno,no3,(p+2),(p+3))])%*%beta_hat-what*(xxs$m[q]-2*xxs$z[q]+1))
                     ,exp(-as.matrix(xxs[q,-c(yno,no3,(p+2),(p+3))])%*%beta_hat +what*(xxs$m[q]-2*xxs$z[q]-1)) ) 
    }
    
    
    xxj           <- cbind(xxs,e)
    xxj$prob      <- 1/(1+e)
    
    no4 <- which( colnames(xxj)=="z" )
    no5 <- which( colnames(xxj)=="m" )
    no6 <- which( colnames(xxj)=="indc" )
    no7 <- which( colnames(xxj)=="e" )
    no8 <- which( colnames(xxj)=="prob" )
    
    multmat       <- xxj[,-c(yno,no4,no5,no6,no7,no8)]*xxj$prob*xxj$indc
    scorevec[i,]  <- apply(multmat,2,sum)
  }
  
  no1 <- which( colnames(xdat)=="id" )
  no2 <- which( colnames(xdat)=="w" )
  
  xdes     <-  xdat[,-c(yno,no1,no2)]
  xdes     <-  as.matrix(xdes)
  
  
  vmat     <-  var(scorevec)
  hmat     <-  t(xdes*(fit$weights))%*%xdes
  Cov_beta <-  n*solve(hmat)%*%vmat%*%solve(hmat)
  
  std_bet  <- sqrt(diag(Cov_beta))
  
  Z        <- beta_hat/std_bet
  pValue   <- 2*pnorm(-abs(Z))  
  
 
  
}

res <- list(beta_hat = beta_hat ,Cov_beta = Cov_beta,  T_stat= T,normal_quantile =quantile )
return(res)

}