HLCR: Highest Likelihood Confidence Region

In ppham27/setsim: Software tool for the visualization of confidence regions in statistical models

Description

This function generates the highest likelihood confidence interval around MLE for a given parameter. It is recommended to use samples from independent sampling methods.

Usage

HLCR(conf, sample, max_loglik, decimal)

Arguments

conf	the desired confidence level
sample	the parameter samples generated by independent sampling method. The last column of this matrix should be the log-likelihood value for each set of parameter
max_loglik	maximum log-likelihood value of the model
decimal	how accurate do you want c to be? default value is 4, which means c is round up to the nearest 0.0001.

Details

This function serves as an add-on to the independent sampling method.

Value

inside	returns parameter samples inside the confidence region
outside	returns parameter samples outside the confidence region
c	the ratio between likelihood of samples and the maximum of the likelihood function in the (1-α) HLCR region for θ.

See Also

independent for more information about independent sampling method.

Examples

# run independent sampling method
library(MASS)
data(data.poisson)
Device <- data.poisson$Device
DataPlan <- data.poisson$DataPlan
y <- data.poisson$y
fit <- glm(y ~ (Device+DataPlan), family=poisson)
out_i <- independent("Mpoisson",y,fit,B=1000)

Region <- HLCR(0.95,out_i$ind[,5:8],out_i$max)
Region$c
Region$inside
Region$outside

# visualize the relationship between c and CR.
conf <- seq(0.01,0.99,by=0.01)
c <- 0
for (i in 1:length(conf)){
  c[i] <- HLCR(conf[i],out_i$ind[,5:8],out_i$max)$c}
plot(c~conf,type="l",xlab="CR",main="Relationship between c and Confidence Region")

##### 3d visualization - Require library(rgl)
# samples inside the confidence region have blue colors
# samples outside the confidence region have green colors
# MLE is red-colored text
library(rgl)
plot3d(x=Region$ins[,2],y=Region$ins[,3],z=Region$ins[,4],xlab="Device",ylab="DataPlan",zlab="log-likelihood",col="blue",xlim=c(min(Region$out[,2]),max(Region$out[,2])),ylim=c(min(Region$out[,3]),max(Region$out[,3])),zlim=c(min(Region$out[,4]),out_i$max))
points3d(x=Region$out[,2],y=Region$out[,3],z=Region$out[,4],col="green")
text3d(x=fit$coef[2],y=fit$coef[3],z=out_i$max,text="MLE",col="red")

ppham27/setsim documentation built on May 25, 2019, 11:25 a.m.