Moment Estimation Methods

1.) Log Density Derivation

Here we are examining the log-logistic (or Fisk) distribution, a popular distribution when dealing with non-negative and rightly skewed data (think survival analysis or income distributions for example).

Using R’s built-in D function, I was able to quickly find this distribution’s relative probability density function (PDF), after already being provided its cumulative distribution function (CDF). These calculations are predicated on the fact that a distribution’s PDF is simply the first derivative of its CDF. Since a derivative is the rate of change of a function at a given point, deriving a CDF across the entirety of its bounds allows for an understanding of how relative probability changes throughout. Below are the Fisk distribution’s PDF and CDF, denoted by f(x) and F(x) respectively.

CDF_function = function(x, alpha, beta){
1/(1+(x/alpha)^(-beta))
}

CDF_expression = expression(
  1/(1+(x/alpha)^(-beta))
  )


# Creating both an expression object and function object since the deriv function I am about to use only takes expression objects
  
  # Will use D function to see the returned PDF in algebraic form
  # Will use deriv function to get an expression that can be used for evaluations

## PDF = first derivative of CDF
PDFD = D(CDF_expression, "x")
#PDFD  # = ((x/alpha)^((-beta) - 1) * ((-beta) * (1/alpha))/(1 + ((x/alpha)^(-beta)))^2)
  # Wrote this down, did some additional algebra before writing up Latex below

PDF_Function = function(x, alpha, beta){
  ifelse(x <= 0, 0,
         (beta/alpha) * (x/alpha)^(beta - 1) / (1 + (x/alpha)^beta)^2)
}
  # Had to include the ifelse component because x will always be positive in the scenario of our data (recovery times) and in accordance with the properties of the Fisk distribution
    # I previously did NOT have that ifelse component and my visual confirmation further along was not effective


PDF_deriv = deriv(CDF_expression, "x")

\[{\Large F(x) = \frac{1}{1 + (x/\alpha)^{-\beta}}} \hspace{0.5cm}\Large , \hspace{0.5cm} X>=0 \]

\[{\Large f(x) = \frac{\beta}{\alpha} * \frac{(x/\alpha)^{\beta-1}}{\alpha(1+(x/\alpha)^\beta)^2}} \hspace{0.5cm}\Large , \hspace{0.5cm}X>=0 \]

2.) Recovery Time Distribution

The dataset we have on hand is composed of 50 different patients’ observed recovery times (in days) following a particular knee surgery procedure. We will operate the analysis of this dataset under the assumption that the underlying population of knee surgery recovery times is distributed in a log-logistic manner. Below are both the summary statistics and a histogram of the dataset’s values.

Days = c(8.23, 12.74, 14.83, 16.61, 18.16, 19.55, 20.80, 21.94, 23.00, 23.98, 24.89, 25.75, 26.56, 
27.34, 28.08, 28.79, 29.48, 30.15, 30.81, 31.45, 32.08, 32.70, 33.31, 33.92, 34.53, 35.13, 
35.73, 36.33, 36.93, 37.53, 38.14, 38.75, 39.37, 40.00, 40.64, 41.29, 41.95, 42.63, 43.33, 
44.05, 44.79, 45.56, 46.36, 47.20, 48.08, 49.02, 50.03, 51.12, 52.32, 53.65)
  # n = 50
Data_Days = data.frame(Days)

sum_stats_table = function(variable){
  table = cbind(
    round(min(variable),2),
    round(quantile(variable, 0.25),2),
    round(median(variable),2),
    round(mean(variable),2),
    round(quantile(variable, 0.75),2),
    round(max(variable),2)
  )
  colnames(table) = c("Min", "Q1", "Med", "Mean", "Q3", "Max")
  rownames(table) = NULL
  return(table)
}
  # I created a table for summary statistics of a variable, can use this moving forward instead of manually putting the summary stats in table form


Summary_Table = sum_stats_table(variable = Days)
kable(Summary_Table,align = "c",
      caption = "<span style='color:##000000;'>Summary Stats of Recovery Times</span>") %>%
  kable_styling(
    bootstrap_options = c("striped", "bordered"),
    full_width = FALSE,
    position = "center")

Summary Stats of Recovery Times

Min	Q1	Med	Mean	Q3	Max
8.23	26.75	34.83	34.19	42.46	53.65

hist(Days, 
     main = "Distribution of Recovery Times",
     xlab = "Recovery Time (Days)",
     ylab = "Count",
)
box(lwd=3)

2a.) Moment Estimation

As stated above, we are operating under the assumption of our data stemming from a Fisk distribution. Therefore, there are two population parameters of interest, \(\alpha\) (rate of scale) and \(\beta\) (rate of shape), which we would like to use our sample data to estimate.

Using moment estimation, general formula tells us that the kth moment (\(m_k\)) of a distribution can be found by summing said dataset’s values (each raised to the k power) and then divided by the sample size. Additionally, and more specifically to the Fisk distribution, we can directly equate the different moments of a distribution to functions of the \(\alpha\) and \(\beta\) parameters. The formulas and calculations of the first two moments (\(m_1\) and \(m_2\)) are below.

## Reference Section 4.2 example 2 from lecture notes

# Mean of log logistic distribution = (alpha * pi)
                                      #      / (beta * sin(alpha/beta))

## Reference formula from Section 3.1
# moment 1 = mean of data # = 34.1922 = (alpha * pi)
                                      #      / (beta * sin(alpha/beta))

# moment 2 = mean of (data squared) = 1288.845

Moment1 = mean(Days)
Moment2 = mean(Days^2)

\[\Large Moment 1 = {\alpha^1*\Gamma(1 + 1/\beta) * \Gamma(1 - 1/\beta)} = \bar{x} \approx 34.1922 \]

\[\Large Moment 2 = {\alpha^2*\Gamma(1 + 2/\beta) * \Gamma(1 - 2/\beta)} = \bar{x^2} \approx 1288.84483 \]

Utilizing our dataset’s first two moments, we then create a system of equations via R’s uniroot and custom function capabilities, to estimate \(\alpha\) and \(\beta\) though numerical methods. Numerical methods are needed since pure algebraic tactics are unable to solve the inequalities for this distribution. Since we are finding estimates of the broader population parameters of a Fisk distributed population, we label these estimates \(\hat{\alpha}\) and \(\hat{\beta}\) to represent \(\alpha\) and \(\beta\) respectively.

## From assignment instructions:
## (alpha^k) * beta function (1 + (k/beta), 1 - (k/beta))
  # Capital B denotes the beta function


Beta_Hat = uniroot(function(Beta){
    # The uniroot function tells us what value of specific variable will lead our function to equal zero or approach zero
    # In this case, we are finding f(beta) = 0 because there is not an algebraic way to find a difference of zero between sample variance and theoretical variance
  
  A = (pi/Beta)/
        sin(pi/Beta)
    # One portion of Fisk Distribution variance formula
  
  B = (2*pi/Beta)/
      sin(2*pi/Beta)
    # Other portion of Fisk Distribution variance formula
  
  Alpha_Beta = Moment1 / A
  Alpha_Beta^2 * (B - A^2) - (Moment2 - Moment1^2)
}, lower = 2.000001, upper = 2000000)$root
  # Set lower bound of search to 2.000001 because, per properties of Fisk Distribution moments, beta must be greater than moment 2 in order for variance to exist
  # upper bound set arbitrarily

Alpha_Hat = Moment1 / ((pi/Beta_Hat)/sin(pi/Beta_Hat))

invisible(
  Alpha_Hat
)

invisible(
  Beta_Hat
)

\[\Large \hat{\alpha} \approx 32.65\] \[\Large \hat{\beta} \approx 6.01\]

To confirm relative accuracy of our estimates, below I overlaid the Fisk distribution’s PDF with given alpha and beta parameters over top our available sample data of patient recovery times. There is moderate overlap between our distribution curve and histogram data. This leads me to believe that there very well could be an underestimate of our \(\alpha\) parameter, as the value of \(\alpha\) traditionally controls the horizontal location of a curve’s peak and decline.

hist(Days, 
     main = "Distribution of Recovery Times",
     xlab = "Recovery Time (Days)",
     ylab = "Probability Density",
     ylim = c(0, 0.05),
     breaks = 10,
     prob = TRUE)


curve(PDF_Function(x , alpha = 32.65431, beta = 6.006248), 
      col = "red",
      lwd = 3,
      add = TRUE)

legend("topright",
       legend = "Fisk (32.65, 6.07)",
       col = "red",
       lwd = 3,
       box.lwd=3)

box(lwd=3)

2b.)

Given the inherent variation that exists in sampling, I chose to execute a bootstrapping sequence. I performed 1,000 iterations of samples (n = 50 each time) for both our parameters estimators, \(\hat{\alpha}\) and \(\hat{\beta}\).

Pictured below are the sampling distributions for both parameters, each of them with a Gaussian kernel density curve layered above and said curve’s respective bandwidth notated in parentheses. Bandwidths were chosen in accordance with Silverman’s standard formula and calculated via a custom-built R function.

# Bootstrap Calculations
set.seed(100)

Alpha_Bootstraps = 0
Beta_Bootstraps = 0
  # Set these variables equal to zero as they will be constantly updated throughout the for loop's iterations

for(i in 1:1000){
  ith_Sample = sample(Days, size = length(Days), replace = TRUE)
    # Always use replacement when bootstrapping
  ith_Moment1 = mean(ith_Sample)
  ith_Moment2 = mean(ith_Sample^2)
  
# Copied function from part 2a: removed comments as process is the same for each bootstrap sample
    # updated variable names for the bootstrap procedure
    Boot_Beta_Hat = uniroot(function(Beta){
    A = (pi/Beta)/
          sin(pi/Beta)
    
    B = (2*pi/Beta)/
        sin(2*pi/Beta)
    
    Alpha_Beta = ith_Moment1 / A
    Alpha_Beta^2 * (B - A^2) - (ith_Moment2 - ith_Moment1^2)
  }, lower = 2.000001, upper = 2000000)$root
  
  Boot_Alpha_Hat = ith_Moment1 / ((pi/Beta_Hat)/sin(pi/Beta_Hat))
  
  Alpha_Bootstraps[i] = Boot_Alpha_Hat
  Beta_Bootstraps[i] = Boot_Beta_Hat
}

## Create custom function for finding Gaussian bw quickly in future
Gaussian_BW = function(Data){
  IQR = quantile(Data, 0.75) - quantile(Data, 0.25)
  Min_Part = min(sd(Data),(IQR/1.34))
  BW = 0.9 * Min_Part * (length(Data)^(-0.2))
  return(BW)
}

invisible(
  Gaussian_BW(Data = Alpha_Bootstraps) # = 0.3183099
)


hist(Alpha_Bootstraps,
     main = "Bootstrap Distribution of Alpha Hat",
     xlab = "Sampling Means",
     ylab = "Probability Density",
     prob = TRUE
     )

# bandwidth calculations below
Gaussian = density(Alpha_Bootstraps, bw = 0.3183099, kernel = "gaussian")
lines(Gaussian,
      lwd = 3,
      col = "blue"
      )

legend("topright",
       legend = "Gaussian (0.32)",
       col = "blue",
       lwd = 3,
       box.lwd=3)

box(lwd=3)

invisible(
  Gaussian_BW(Data = Beta_Bootstraps) # = 0.1303286
)

hist(Beta_Bootstraps,
     main = "Bootstrap Distribution of Beta Hat",
     xlab = "Sampling Means",
     ylab = "Probability Density",
     prob = TRUE,
     ylim = c(0,0.7)
     )
box(lwd=3)

Gaussian = density(Beta_Bootstraps, bw = 0.1303286, kernel = "gaussian")
lines(Gaussian,
      lwd = 3,
      col = "blue"
      )

legend("topright",
       legend = "Gaussian (0.13)",
       col = "blue",
       lwd = 3,
       box.lwd=3)

What the bootstrap distribution for \(\hat{\alpha}\) and \(\hat{\beta}\) display is a far greater resemblance with their Gaussian PDF estimations than what was seen above, comparing our data to the curve of a standard Fisk distribution with said values for \(\alpha\) and \(\beta\).

Two possible suggestions that this relevation brings are: perhaps the population from which these sampled surgery recovery times came from is not truly as Fisk distributed as we initially believed, or as stated previously there could be an underestimate of the \(\alpha\) parameter.