r - multilevel zinb.txt

R 多变量分析工具

	{
		if(model != "rlg" & model != "zinb")
	{
		glm.poi <- wreml.poi(y,wmm,x.p,rv,beta,va,sig2)
		beta <- glm.poi$beta
		va <- glm.poi$va
	}
		else
		beta <- coef(glm(y ~ x.p, family = poisson(link =log),
					       weights = wmm, na.action = na.omit, control = ct0))
	}
	else 
	beta <- coef(glm(y ~ 1, family = poisson(link =log), 
					       weights = wmm, na.action = na.omit, control = ct0))
	
	if(max(abs(ZK-ZK1))<1e-3) {flag <- 1;break;}
   	ZK1 <- ZK	
cat('\n',iter,k,alfa,'\n')
} # end of inner loop

if(is.null(x.l)) pai <- mean(ZK1)
else pai  <-  theta/(1 + theta)
	
	if(model=="zinbmix")
	{  
		hz <- hznb(y,X,G,rv,pai,mu,th,sig2,sigu)
		vd <- hz$dd
		nsigu  <- as.numeric(t(yu)%*%yu+sum(vd[pa+pb+1:m]))/m
		nsigv2  <- as.numeric(t(va)%*%va+sum(vd[pa+pb+m+1:m]))/m
		if(nsigv2<=0) nsigv2 <- 1
		if(nsigu<=0) nsigu <- 1
	#	cat("model=zinbmix")	
	}
	if(model=="rnb")
	{	vd <- hznb(y,X,G,rv,pai,mu,th,sig2)
		nsigv2  <- as.numeric(t(va)%*%va+sum(vd[(pa+pb+1):(pa+pb+m)]))/m
		if(nsigv2<=0)
		nsigv2 <- 1;
		sigu <- nsigu <- 1;
		cat("model=rnb")
    }  
	if(model=="rlg")         
	{	vd <- hznb(y,X,G,rv,pai,mu,th,,sigu)
		nsigu  <- as.numeric(t(yu)%*%yu+sum(vd[(pa+pb+1):(pa+pb+m)]))/m
		if(nsigu<=0) nsigu <- 1;
		sig2 <- nsigv2 <- 1;
		cat("model=rlg")
	    } 
	if(model=="zinb") 
	{
	vd <- hznb(y,X,G,rv,pai,mu,th)	
	flag <- 1;cat("model=zinb")
	break
	}          
if((abs(nsigv2-sig2)<1e-4)&(abs(nsigu-sigu)<1e-4)) {flag <- 1;break}
	sig2 <- nsigv2	
	sigu <- nsigu
cat('\n',ie,sig2,sigu,th,'\n')
} #end of loop

names(beta) <- c("Intercept", dimnames(x.p)[[2]])  # NOTE: problem here?
names(alfa) <- c("Intercept", dimnames(x.l)[[2]])
std <- sqrt(vd)

if(model=="zinbmix")
{
	risku <- as.vector(yu)#/std[pa+pb+1:m])
	riskv <- as.vector(va)#/std[pa+pb+m+1:m])
	names(riskv) <- names(risku) <- names(table(random))
}
if(model=="rnb")
{
riskv <- as.vector(va/std[pa+pb+1:m])
names(riskv) <- names(table(random))
#cat("OK2")
}
if(model=="rlg")
{
	risku <- as.vector(yu/std[pa+pb+1:m])
	names(risku) <- names(table(random))
}

stda <- std[1:pa]
stdb <- std[pa+1:pb]

eta <- cbind(alfa,stda,alfa/stda,2*(1-pnorm(abs(alfa/stda))))
dimnames(eta) <- list(names(alfa),c("Estimate","SD","t-value","p-value"))
etb <- cbind(beta,stdb,beta/stdb,2*(1-pnorm(abs(beta/stdb))))
dimnames(etb) <- list(names(beta),c("Estimate","SD","t-value","p-value"))
eta <- round(eta,3)
etb <- round(etb,3)

obj.call <- match.call()
if(model=="zinbmix")
{
	result <- list(call=obj.call, th=th, pai=pai, mu=mu, beta=etb, alfa=eta)
	result$sigv <- round(sig2,4)
	result$riskv <- round(riskv,4)
	result$sigu <- round(sigu,4)
	result$risku <- round(risku,4)
	result$df <- hz$df
	#cat("OKmix")

}
if(model=="rnb")
{
result <- list(call=obj.call, th=th, beta=etb, pai=pai, mu=mu, alfa=eta, 
                sigv=sig2, riskv=round(riskv,4))
cat("OKrnb")
}
if(model=="rlg")
{	result <- list(call=obj.call, th=th, pai=pai, mu=mu, beta=etb, alfa=eta)
	result$sigu <- round(sigu,4)
	result$risku <- round(risku,4)
	cat("OKrlg")

}
if(model=="zinb")
{result <- list(call=obj.call, th=th, beta=etb, pai=pai, mu=mu, alfa=eta)
cat("OKzinb")
}
if(flag) result
else stop("error:not reach the convergence")
}

#-----------------
# main function
#-----------------
# function sumzinbmix()
# y------count
# x.nb-----covariate matrix for mean
# x.l ----covariate matrix for logistic part
# random, r.nb,r.l are flags for whether including random effects in the models.

sumzinbmix <- function(y, x.nb=NULL, x.l=NULL, random=NULL, r.nb=F, r.l=F)
{
	if(!is.null(random))
	{Group <- c(as.factor(random)) # DAVE: codes --> c
		n <- length(y)
		m <- max(Group)
		z <- matrix(0, ncol = m, nrow = n)
		for(i in 1:m) 
		z[, i] <- ifelse(Group == i, 1,0)
	}
	else z <- diag(10)
if(r.nb&r.l)
	obj <- zinbmix(y, x.nb, z, random, x.l, model="zinbmix")	# random effects in both parts
if(r.nb&(!r.l))
	obj <- zinbmix(y, x.nb, z, random, x.l, model="rnb")      # random in nb part	
if((!r.nb)&r.l)
	obj <- zinbmix(y, x.nb, z, random, x.l, model="rlg")      # random in logistic part	
if((!r.nb)&(!r.l))
	obj <- zinbmix(y, x.nb, z, random, x.l, model="zinb")     # no random 	

#read the object
	pai <- obj$pai
	th <- obj$th
	beta <- obj$beta
	alfa <- obj$alfa #parameters in logistic part
	mu <- obj$mu #mean of y
	
#*****************************************************************

#*****************************************************************
#--------------calculate the sumarry results  -------------------	 #frequency

k <- max(y) ; n <- length(y)
fr.z <- fr.ob <- 0:k
fr.z[1] <- sum(pai+(1-pai)*(th/(th+mu))^th)
for(jj in 1:k){
fr.z[jj+1] <- sum((gamma(jj+th)/gamma(jj+1)/gamma(th)*(th/(mu+th))^th*(mu/(mu+th))^jj)*(1-pai))
}
for(jj in 0:k){fr.ob[jj+1] <- sum(ifelse(y==jj,1,0))}
fr.ob <- round(fr.ob[1:(k+1)],3);
fr.z <- round(fr.z[1:(k+1)],3)
devia <<- cbind(0:k,(fr.z-fr.ob)/n)

fr.with <- matrix(c((0:k),fr.ob,fr.z),ncol=3)
dimnames(fr.with) <- list(NULL,c("count","observed","expected"))


fr.ob <- c(fr.ob[fr.z>5],n-sum(fr.ob[fr.z>5]))
fr.z <- c(fr.z[fr.z>5],n-sum(fr.z[fr.z>5]))
chi.z <- sum((fr.ob-fr.z)^2/fr.z)
pv <- 1-pchisq(chi.z,length(fr.z)-4)

#std error of alpha

ksi <- pai/(1-pai)
u <- th/(th+mu)
yzero <- ifelse(y > 0, 0, 1)
#____________________________________
f0 <- table(y[y > 0])
f <- rep(0, k)
f[as.numeric(names(f0))] <- f0
tot <- sum(f0)
f <- tot + f - cumsum(f)
#------------------------------------	

# first deriv of A(th)
i <- sum(f/(th+1:k-1))
ii <- sum(f/(th + 1:k-1)^2)

B <- log(u)+(1-u)
B1 <- (1-u)^2/th
ep1 <- u^th
ep2 <- (ksi+ep1)^2

w33 <- sum(yzero*ksi*(B1*(ksi+ep1)-B^2*ep1)/ep2-B1-(1-yzero)*y*(u/th)^2)+ii
sdc <- sqrt(1/w33)/th^2
etc <- cbind(1/th,sdc,1/sdc/th,2*(1-pnorm(abs(1/sdc/th))))
etc <- round(etc,4)
dimnames(etc) <- list("alpha",c("Estimate","SD","t-value","p-value"))

# log-likehood

loglik <- (sum(yzero*log(ksi+ep1)-log(1+ksi)+(1-yzero)*(lgamma(y+th)-lgamma(y+1)-lgamma(th)+log(ep1)+y*log(1-u))))

#pearson residuals

mu.y <- (1-pai)*mu
std.y <- sqrt(mu.y*(1+(1/th+pai)*mu))
r.p <- (y-mu.y)/std.y 
PearChi <- round(sum(r.p^2),3)   #Pearson residuals
r.p <- round(r.p,4)   #residuals

#Deviance residual

lnb <- (lgamma(y+th)-lgamma(y+1)-lgamma(th)+th*log(th/(y+th))+log((y/(th+y))^y))

lzinb <- (yzero*log(ksi+ep1)-log(1+ksi)+(1-yzero)*(lgamma(y+th)-lgamma(y+1)-lgamma(th)+log(ep1)+y*log(1-u)))

d <- 2*(lnb-lzinb)
Dev <- round(sum(d),3)# Deviance;
r.d <- round(sign(y-mu.y)*sqrt(d),4) #Deviance residuals

residm <<- cbind(y, mu.y, r.p, r.d, lnb, lzinb)    # output to a data sheet

#-----------------------------------------------------
#******************************************************
#******************************************************
# setTextOutputRouting("Report","Default") # DAVE: escape text routing
cat("\n","Fit of ZINB mixed model",'\n')	
cat('\n',"inflated part",'\n')
print(obj$alfa)
cat('\n','_____________________________________','\n\n')

cat('\n\n',"negative binomial part",'\n')
print(obj$beta)
cat('\n','_____________________________________','\n\n')
print(etc)

if(r.l)
{
cat('\n',"sigma^2 of random effect in inflate :", round(obj$sigu,5),'\n')
cat('\n',"random effects:\n")
risku <<- as.matrix(obj$risku)
print(obj$risku)
cat('\n','_____________________________________','\n\n')

}
if(r.nb)
{
cat('\n',"sigma^2 of random effect in NB :", round(obj$sigv,5),'\n')
cat('\n',"random effects:\n")
riskv <<- as.matrix(obj$riskv)   #
print(obj$riskv)
cat('\n','_____________________________________','\n\n')
}

cat('\n',"--------------------------------------",'\n\n')
obj$count.tab <- fr.with      # Dave: add table of counts to output object
print(fr.with)
cat('\n','_____________________________________','\n\n')
cat('\n\n',"Chisquare test statistics: ",round(chi.z,3))
cat('\n',"loglikelihood              : ",round(loglik,3))
cat('\n',"Pearson residuals          : ",round(PearChi,3))
obj$chi.sq <- chi.z           # Dave: save indices to model object
obj$loglik <- loglik
# cat('\n',"Deviance                   : ",round(Dev, 3))
# cat('\n', round(obj$df,3)) # Dave: error msg ?

# setTextOutputRouting("Default","Default") # Dave: escape text-outputting
# cat("\n")
# return("End of program.  See results in the report window !") # Dave: escape
#
### Dave: output whole model object
return(obj)
}

# example
# sumzinbmix(y,x.nb=NULL,,x.l=NULL,random=NULL,r.nb=F,r.l=F)

######################### Coefficient Extraction Function ######################

zinb.coef <- function(obj, dig=2){
  lr.coef <- data.frame(obj$alfa)
  nb.coef <- data.frame(obj$beta)
  lr.coef$OR <- exp(lr.coef[,"Estimate"])
  lr.coef$Lower <- exp(lr.coef[,"Estimate"] - 1.96*lr.coef[,"SD"])
  lr.coef$Upper <- exp(lr.coef[,"Estimate"] + 1.96*lr.coef[,"SD"])
  nb.coef$exp.b <- exp(nb.coef[,"Estimate"])
  nb.coef$Lower <- exp(nb.coef[,"Estimate"] - 1.96*nb.coef[,"SD"])
  nb.coef$Upper <- exp(nb.coef[,"Estimate"] + 1.96*nb.coef[,"SD"])
  cat("\n", "Logistic Model", "\n")
  print(round(lr.coef, dig))
  cat("\n", "Negative Binomial Model", "\n")
  print(round(nb.coef, dig))
  }