diff --git a/R/predict.hmsc.R b/R/predict.hmsc.R
index 6eff434..9841bd3 100644
--- a/R/predict.hmsc.R
+++ b/R/predict.hmsc.R
@@ -496,22 +496,103 @@ predict.hmsc<-function(object, newdata, type = c("response","link"), conditional
 		### Apply inverse link function
 		if(type=="response"){
 			if(any(class(object)=="probit")){
-				result<-pnorm(apply(res,1:2, mean))
+				#result<-pnorm(apply(res,1:2, mean))
+mu <- apply(res, 1:2, mean)
+q025 <- apply(res, 1:2, function(x) {
+  pnorm(quantile(x, .025))
+})
+q05 <- apply(res, 1:2, function(x) {
+  pnorm(quantile(x, .05))
+})
+q10 <- apply(res, 1:2, function(x) {
+  pnorm(quantile(x, .1))
+})
+q90 <- apply(res, 1:2, function(x) {
+  pnorm(quantile(x, .9))
+})
+q95 <- apply(res, 1:2, function(x) {
+  pnorm(quantile(x, .95))
+})
+q975 <- apply(res, 1:2, function(x) {
+  pnorm(quantile(x, .975))
+})
+result <- list(mu = mu, q025 = q025, q05=q05, q10 = q10, q90 = q90, q95 = q95, q975 = q975)
 			}
 
 			if(any(class(object)=="gaussian")){
-				result<-apply(res,1:2, mean)
+				#result<-apply(res,1:2, mean)
+mu <- apply(res, 1:2, mean)
+q025 <- apply(res, 1:2, function(x) {
+  quantile(x, .025)
+})
+q05 <- apply(res, 1:2, function(x) {
+  quantile(x, .05)
+})
+q10 <- apply(res, 1:2, function(x) {
+  quantile(x, .1)
+})
+q90 <- apply(res, 1:2, function(x) {
+  quantile(x, .9)
+})
+q95 <- apply(res, 1:2, function(x) {
+  quantile(x, .95)
+})
+q975 <- apply(res, 1:2, function(x) {
+  quantile(x, .975)
+})
+result <- list(mu = mu, q025 = q025, q05=q05, q10 = q10, q90 = q90, q95 = q95, q975 = q975)
 			}
 
 			if(any(class(object)=="poisson" | any(class(object)=="overPoisson"))){
-				result<-exp(apply(res,1:2, mean))
+				#result<-exp(apply(res,1:2, mean))
+mu <- apply(res, 1:2, mean)
+q025 <- apply(res, 1:2, function(x) {
+  exp(quantile(x, .025))
+})
+q05 <- apply(res, 1:2, function(x) {
+  exp(quantile(x, .05))
+})
+q10 <- apply(res, 1:2, function(x) {
+  exp(quantile(x, .1))
+})
+q90 <- apply(res, 1:2, function(x) {
+  exp(quantile(x, .9))
+})
+q95 <- apply(res, 1:2, function(x) {
+  exp(quantile(x, .95))
+})
+q975 <- apply(res, 1:2, function(x) {
+  exp(quantile(x, .975))
+})
+result <- list(mu = mu, q025 = q025, q05=q05, q10 = q10, q90 = q90, q95 = q95, q975 = q975)
 			}
 		}
 
 		if(type=="link"){
-				result<-apply(res,1:2, mean)
+				#result<-apply(res,1:2, mean)
+mu <- apply(res, 1:2, mean)
+q025 <- apply(res, 1:2, function(x) {
+  quantile(x, .025)
+})
+q05 <- apply(res, 1:2, function(x) {
+  quantile(x, .05)
+})
+q10 <- apply(res, 1:2, function(x) {
+  quantile(x, .1)
+})
+q90 <- apply(res, 1:2, function(x) {
+  quantile(x, .9)
+})
+q95 <- apply(res, 1:2, function(x) {
+  quantile(x, .95)
+})
+q975 <- apply(res, 1:2, function(x) {
+  quantile(x, .975)
+})
+result <- list(mu = mu, q025 = q025, q05=q05, q10 = q10, q90 = q90, q95 = q95, q975 = q975)
 		}
-		colnames(result)<-colnames(Y)
+		#colnames(result)<-colnames(Y)
+	  ## the colnames match without this
 	}else{
 		### Apply inverse link function
 		if(type=="response"){