From bb56394865a3cebfed309b4829ca772c1a1972a0 Mon Sep 17 00:00:00 2001 From: sastoudt Date: Tue, 9 Apr 2019 16:18:54 -0700 Subject: [PATCH 1/3] return various quantiles for uncertainty quantification in prediction --- R/predict.hmsc.R | 47 ++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 42 insertions(+), 5 deletions(-) diff --git a/R/predict.hmsc.R b/R/predict.hmsc.R index 6eff434..2960aad 100644 --- a/R/predict.hmsc.R +++ b/R/predict.hmsc.R @@ -496,22 +496,59 @@ 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)) + resultInt= apply(res,1:2, function(x){list(mu=mean(x),q05=quantile(x,.05),q10=quantile(x,.1),q90=quantile(x,.9),q95=quantile(x,.95))}) + 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,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) + resultInt= apply(res,1:2, function(x){list(mu=mean(x),q05=quantile(x,.05),q10=quantile(x,.1),q90=quantile(x,.9),q95=quantile(x,.95))}) + 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,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)) + resultInt= apply(res,1:2, function(x){list(mu=mean(x),q05=quantile(x,.05),q10=quantile(x,.1),q90=quantile(x,.9),q95=quantile(x,.95))}) + 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,q10=q10,q90=q90,q95=q95,q975=q975) ## fix here } } if(type=="link"){ - result<-apply(res,1:2, mean) + #result<-apply(res,1:2, mean) + resultInt= apply(res,1:2, function(x){list(mu=mean(x),q05=quantile(x,.05),q10=quantile(x,.1),q90=quantile(x,.9),q95=quantile(x,.95))}) + 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,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"){ From dd3ae51b160d5dd84b3315c812f622273351ccc4 Mon Sep 17 00:00:00 2001 From: sastoudt Date: Tue, 9 Apr 2019 16:23:42 -0700 Subject: [PATCH 2/3] clean up --- R/predict.hmsc.R | 128 ++++++++++++++++++++++++++++++++++------------- 1 file changed, 92 insertions(+), 36 deletions(-) diff --git a/R/predict.hmsc.R b/R/predict.hmsc.R index 2960aad..136f3f7 100644 --- a/R/predict.hmsc.R +++ b/R/predict.hmsc.R @@ -497,55 +497,111 @@ predict.hmsc<-function(object, newdata, type = c("response","link"), conditional if(type=="response"){ if(any(class(object)=="probit")){ #result<-pnorm(apply(res,1:2, mean)) - resultInt= apply(res,1:2, function(x){list(mu=mean(x),q05=quantile(x,.05),q10=quantile(x,.1),q90=quantile(x,.9),q95=quantile(x,.95))}) - 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,q10=q10,q90=q90,q95=q95,q975=q975) + resultInt <- apply(res, 1:2, function(x) { + list(mu = mean(x), q05 = quantile(x, .05), q10 = quantile(x, .1), q90 = quantile(x, .9), q95 = quantile(x, .95)) +}) +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, q10 = q10, q90 = q90, q95 = q95, q975 = q975) } if(any(class(object)=="gaussian")){ #result<-apply(res,1:2, mean) - resultInt= apply(res,1:2, function(x){list(mu=mean(x),q05=quantile(x,.05),q10=quantile(x,.1),q90=quantile(x,.9),q95=quantile(x,.95))}) - 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,q10=q10,q90=q90,q95=q95,q975=q975) + resultInt <- apply(res, 1:2, function(x) { + list(mu = mean(x), q05 = quantile(x, .05), q10 = quantile(x, .1), q90 = quantile(x, .9), q95 = quantile(x, .95)) +}) +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, q10 = q10, q90 = q90, q95 = q95, q975 = q975) } if(any(class(object)=="poisson" | any(class(object)=="overPoisson"))){ #result<-exp(apply(res,1:2, mean)) - resultInt= apply(res,1:2, function(x){list(mu=mean(x),q05=quantile(x,.05),q10=quantile(x,.1),q90=quantile(x,.9),q95=quantile(x,.95))}) - 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,q10=q10,q90=q90,q95=q95,q975=q975) ## fix here + resultInt <- apply(res, 1:2, function(x) { + list(mu = mean(x), q05 = quantile(x, .05), q10 = quantile(x, .1), q90 = quantile(x, .9), q95 = quantile(x, .95)) +}) +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, q10 = q10, q90 = q90, q95 = q95, q975 = q975) } } if(type=="link"){ #result<-apply(res,1:2, mean) - resultInt= apply(res,1:2, function(x){list(mu=mean(x),q05=quantile(x,.05),q10=quantile(x,.1),q90=quantile(x,.9),q95=quantile(x,.95))}) - 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,q10=q10,q90=q90,q95=q95,q975=q975) + resultInt <- apply(res, 1:2, function(x) { + list(mu = mean(x), q05 = quantile(x, .05), q10 = quantile(x, .1), q90 = quantile(x, .9), q95 = quantile(x, .95)) +}) +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, q10 = q10, q90 = q90, q95 = q95, q975 = q975) } #colnames(result)<-colnames(Y) ## the colnames match without this From f17937eb088637c2e6e18f0079294b1fb28280fc Mon Sep 17 00:00:00 2001 From: sastoudt Date: Tue, 9 Apr 2019 16:37:43 -0700 Subject: [PATCH 3/3] fix redundancy and missing return quantile --- R/predict.hmsc.R | 20 ++++---------------- 1 file changed, 4 insertions(+), 16 deletions(-) diff --git a/R/predict.hmsc.R b/R/predict.hmsc.R index 136f3f7..9841bd3 100644 --- a/R/predict.hmsc.R +++ b/R/predict.hmsc.R @@ -497,9 +497,6 @@ predict.hmsc<-function(object, newdata, type = c("response","link"), conditional if(type=="response"){ if(any(class(object)=="probit")){ #result<-pnorm(apply(res,1:2, mean)) - resultInt <- apply(res, 1:2, function(x) { - list(mu = mean(x), q05 = quantile(x, .05), q10 = quantile(x, .1), q90 = quantile(x, .9), q95 = quantile(x, .95)) -}) mu <- apply(res, 1:2, mean) q025 <- apply(res, 1:2, function(x) { pnorm(quantile(x, .025)) @@ -519,14 +516,11 @@ q95 <- apply(res, 1:2, function(x) { q975 <- apply(res, 1:2, function(x) { pnorm(quantile(x, .975)) }) -result <- list(mu = mu, q025 = q025, q10 = q10, q90 = q90, q95 = q95, q975 = q975) +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) - resultInt <- apply(res, 1:2, function(x) { - list(mu = mean(x), q05 = quantile(x, .05), q10 = quantile(x, .1), q90 = quantile(x, .9), q95 = quantile(x, .95)) -}) mu <- apply(res, 1:2, mean) q025 <- apply(res, 1:2, function(x) { quantile(x, .025) @@ -546,14 +540,11 @@ q95 <- apply(res, 1:2, function(x) { q975 <- apply(res, 1:2, function(x) { quantile(x, .975) }) -result <- list(mu = mu, q025 = q025, q10 = q10, q90 = q90, q95 = q95, q975 = q975) +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)) - resultInt <- apply(res, 1:2, function(x) { - list(mu = mean(x), q05 = quantile(x, .05), q10 = quantile(x, .1), q90 = quantile(x, .9), q95 = quantile(x, .95)) -}) mu <- apply(res, 1:2, mean) q025 <- apply(res, 1:2, function(x) { exp(quantile(x, .025)) @@ -573,15 +564,12 @@ q95 <- apply(res, 1:2, function(x) { q975 <- apply(res, 1:2, function(x) { exp(quantile(x, .975)) }) -result <- list(mu = mu, q025 = q025, q10 = q10, q90 = q90, q95 = q95, q975 = q975) +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) - resultInt <- apply(res, 1:2, function(x) { - list(mu = mean(x), q05 = quantile(x, .05), q10 = quantile(x, .1), q90 = quantile(x, .9), q95 = quantile(x, .95)) -}) mu <- apply(res, 1:2, mean) q025 <- apply(res, 1:2, function(x) { quantile(x, .025) @@ -601,7 +589,7 @@ q95 <- apply(res, 1:2, function(x) { q975 <- apply(res, 1:2, function(x) { quantile(x, .975) }) -result <- list(mu = mu, q025 = q025, q10 = q10, q90 = q90, q95 = q95, q975 = q975) +result <- list(mu = mu, q025 = q025, q05=q05, q10 = q10, q90 = q90, q95 = q95, q975 = q975) } #colnames(result)<-colnames(Y) ## the colnames match without this