# R/BIest.R In HRQoL: Health Related Quality of Life Analysis

```BIest <- function(y,m,disp=FALSE,method="MM"){

if (sum(as.integer(y)==y)==length(y)){
} else {
stop("y must be integer")
}

if (sum(m==as.integer(m))==length(m)){
} else {
stop("m must be integer")
}

if ((length(m)==1) | (length(m)==length(y))){
} else{
stop("m must be a number, or a vector of the length of y")
}

if (max(y-m)>0 | min(y) < 0){
stop("y must be bounded between 0 and m")
}

if (min(m)<=0){
stop("m must be positive")
}

if ((disp==FALSE) | (disp==TRUE)){
} else {
stop("disp is logical")
}

if ((method=="MM") | (method=="MLE")){
} else{
stop("method option must be MM or MLE")
}

#Number of observations
n <- length(y)
#Balanced data
if (length(m)==1){
m <- rep(m,n)
balanced <- "yes"
} else{
if(sum(m==m)==length(m)){
balanced <- "yes"
} else {
balanced <- "no"
}
}

# The estimating equation for p is the same with both methodologies!
p <- sum(y)/sum(m)
# Ones we have computed the mle, we replace it in the Fisher information formula
I <- sum(m)/(p*(1-p))

if (disp){
if (n==1){
stop("Warning: Cannot calculate overdipersion with one observation")
}
a <- y==y
if (sum(a)==n){
stop("Warning: Cannot calculate overdispersion if all the observations are equal")
}

if (method=="MLE"){
#Estimate the dispersion parameter
y0. <- which(y==0)
if (length(y0.)==0){
phi0 <- 0
} else{
y0 <- y[y0.]
m0 <- m[y0.]
phi0 <- sum(m0*log(1/(1-p)))
}

ym. <- which(y==m)
if (length(ym.)==0){
phim <- 0
} else {
ym <- y[ym.]
mm <- m[ym.]
phim <- sum(mm*log(1/p))
}

y0m. <- which((y==0) | (y==m))
if (length(y0m.)==0){
y0m <- y
m0m <- m
phi0m <- sum(y0m*log((y0m/m0m)/p)+(m0m-y0m)*log((1-y0m/m0m)/(1-p)))
} else {
y0m <- y[-y0m.]
m0m <- m[-y0m.]
phi0m <- sum(y0m*log((y0m/m0m)/p)+(m0m-y0m)*log((1-y0m/m0m)/(1-p)))
}

phi <- (2/(n-1))*(phi0+phi0m+phim)
} else {
phi <- sum((y-m*p)^2/(m*p*(1-p)))/(n-1)
}
} else {
phi <- 1
}

#Stand. Errors of p
pVar <- phi/I
pSE <- sqrt(pVar)
pIC <- c(p-1.96*pSE,p+1.96*pSE)

out <- list(p=p, pVar=pVar, pIC.low=pIC,pIC.up=pIC,
phi=phi,
m=m,balanced=balanced,
method=method)

class(out) <- "BIest"

out\$call <- match.call()

out
}

print.BIest <- function(x,...){

cat("The probability parameter\n")
cat("Estimation:",x\$p,"\n")
cat("Standard deviation:",sqrt(x\$pVar),"\n")
cat("95% confidence interval:", c(x\$pIC.low,x\$pIC.up),"\n")
if (x\$phi!=1){
cat("\nThe dispersion parameter\n")
cat("Estimation:",x\$phi,"\n")
}
if (x\$balanced=="yes"){
cat("\nBalanced data, maximum score in the trials:",x\$m)
} else {
cat("\nNo balanced data.")
}
cat("\nEstimation approach:",x\$method)
}
```

## Try the HRQoL package in your browser

Any scripts or data that you put into this service are public.

HRQoL documentation built on May 2, 2019, 5:42 a.m.