解釋與分析

1. QQ Plot 的意義：
   • QQ 圖用於檢查數據是否符合正態分佈。
   • 理論正態分佈的分位數由紅線表示，如果數據點（黑色圓點）與紅線接近或重合，則數據分佈接近正態分佈。
2. 數據分佈特徵：
   • 中部區域：
     • 大部分 TPY 數據點位於紅線附近，表明 TPY 數據在中間區域接近正態分佈。
   • 尾部區域：
     • 圖中的左下角和右上角的數據點偏離紅線，表明 TPY 數據在分佈的尾部存在偏差。
     • 右上方的數據點偏高，可能暗示數據分佈右偏（正偏態）。
3. 整體評估：
   • TPY 數據整體上符合正態分佈假設，但尾部（特別是右尾）出現了異常值或偏離現象。

結論

與原始數據 QQ 圖的比較：

• 相較於未取對數時的 TPY 數據，取 ln 後的 LOGTPY 數據在尾部的偏差大幅降低。
• 整體分佈更加對稱且接近正態分佈。
• 取 ln 後的 LOGTPY 數據呈現良好的正態性，大部分數據點均接近紅色理論線。
• 尾部的輕微偏離可以被接受，特別是對多數統計分析而言，取對數後的數據更適合用於後續建模與檢定。

目的是為了改善數據特性，包括減少偏態、減弱異常值影響、轉換非線性關係，以及提高模型的穩定性和解釋力。針對 TPY 數據，取對數是合理且有效的步驟，為後續分析提供了更好的數據基礎。

模型摘要：

1. 係數：
   • 截距：-0.65950，表示當 NUMBED = 0 時，TPY 的預測值。
   • NUMBED 的係數：0.92142，表示 NUMBED 每增加 1 單位，TPY 平均增加 0.92142 單位。
2. 決定係數 (R²)：
   • 模型的決定係數為 R² = 0.9586，這意味著 NUMBED 可以解釋約 95.86% 的TPY方差，表明兩者有一定程度的線性關係。
3. t 統計量：
   • NUMBED 的係數的 t 統計量為 91.346，對應 p 值 < 2e-16。
   • 表明 NUMBED 與 TPY 的關係在統計上顯著。

NUMBED 是 TPY 的一個非常強有力的預測因子，模型能夠很好地解釋 TPY 的變化。

結論

我們比較了兩個線性回歸模型，分別使用 NUMBED 和 SQRFOOT 作為自變數來預測 TPY，並根據模型的決定係數 (R²) 來評估擬合度：

1. 使用 NUMBED 的模型：
   • 決定係數 R² = 0.9586，這意味著 NUMBED 解釋了約 95.86% 的 TPY 方差，顯示出該模型具有較高的擬合度。
2. 使用 SQRFOOT 的模型：
   • 決定係數 R² = 0.6797，這意味著 SQRFOOT 解釋了約 67.97% 的 TPY 方差，顯示出該模型的擬合度相對較低。

根據 R^2 的比較結果，使用 NUMBED 的模型擬合效果更好，能夠解釋更多的 TPY 方差，因此該模型的預測能力較強。

結論：

使用 LOGNUMBED 預測 LOGTPY 的線性回歸模型具有非常高的擬合度，能夠解釋約 94.83% 的變異，並且 LOGNUMBED 對 LOGTPY 的影響在統計上顯著。這表明 LOGNUMBED 是 LOGTPY 的強有力預測因子。

結論：

使用 LOGSQRFOOT 預測 LOGTPY 的回歸模型表現出中等的擬合度，能解釋約 67.65% 的變異，並且 LOGSQRFOOT 對 LOGTPY 的影響在統計上顯著。儘管模型的擬合度不如 LOGNUMBED 模型那麼高，但仍能有效預測 LOGTPY。

線性回歸模型結論

使用 2001 年的數據擬合了以下的線性回歸模型，並檢視其結果。回歸模型的詳細分析如下：

回歸係數

• 截距項：
  • 估計值：-1.034922
  • p 值：0.227
    
  • 結論：截距項在統計上不顯著。
• NUMBED 的回歸係數：
  • 估計值：0.932384
  • p 值：< 2e-16
    
  • 結論：NUMBED 對 TPY 的影響在統計上非常顯著。

擬合度

• 決定係數 \(R^2\)：
  • \(R^2 = 0.9762\)：模型的決定係數表明，NUMBED 解釋了約 97.62% 的 TPY 方差，顯示出模型具有極高的擬合度。
• 調整後的 \(R^2\)：
  • 調整後的 \(R^2 = 0.9762\)：進一步證實了模型的高擬合度。

總結

基於 2001 年的數據，NUMBED 對 TPY 具有顯著的預測能力，並且模型具有非常高的擬合度 （\(R^2 = 0.9762\)）。這表明該模型在 2001 年的數據中表現穩定，且 NUMBED 在預測 TPY 方面的影響未有顯著變化。

Create basic summary statistics for each variable.Summarize the relationship through a correlation statistic and a scatter plot.

1. 基本摘要統計量

library(readr)
data <- read_csv("C:/Users/霍致光/Desktop/大學/應用機率模型/應用期末報告/WiscNursingHome.csv")
library(dplyr)

filtered_data <- data %>% 
  filter(CRYEAR == 2001) %>% 
  select(TPY, NUMBED, SQRFOOT)

filtered_data$LOGTPY <- log(filtered_data$TPY)
filtered_data$LOGNUMBED <- log(filtered_data$NUMBED)

data2 <- filtered_data %>%
select(LOGNUMBED, LOGTPY)

# 1. 摘要統計量
summary(data2)

##    LOGNUMBED         LOGTPY     
##  Min.   :2.890   Min.   :2.511  
##  1st Qu.:4.094   1st Qu.:4.041  
##  Median :4.500   Median :4.396  
##  Mean   :4.457   Mean   :4.368  
##  3rd Qu.:4.779   3rd Qu.:4.700  
##  Max.   :6.125   Max.   :6.088

• 使用 summary() 函數查看了 LOGNUMBED 和 LOGTPY 的基本統計量，包括最小值、第一四分位數、中位數、均值、第三四分位數和最大值。這些摘要統計量有助於了解這兩個變量的分佈情況。

2. 相關係數分析

# 2. 相關係數
correlation <- cor(filtered_data$LOGTPY, filtered_data$LOGNUMBED)
print(correlation)

## [1] 0.9830461

• 計算得到的皮爾森相關係數為 r = 0.9830461

3. 散點圖分析

# 3. 繪製散佈圖
plot(filtered_data$LOGNUMBED, filtered_data$LOGTPY,
     xlab = "LOGNUMBED (床位數的對數)",
     ylab = "LOGTPY (患者年數的對數)",
     main = "LOGNUMBED 與 LOGTPY 的散佈圖",
     pch = 19, col = "blue")

# 添加回歸線
abline(lm(filtered_data$LOGTPY ~ filtered_data$LOGNUMBED), col = "red")

• 散點圖 進一步視覺化了兩者之間的關係，從圖中我們可以清晰地看到，大多數點都大致位於回歸線上，顯示出兩者之間的強線性關係。

Cite the basic summary statistics,include the coefficient of determination, the regression coefficient for LOGNUMBED, and the corresponding t -statistic.

# 擬合線性回歸模型
model <- lm(LOGTPY ~ LOGNUMBED, data = filtered_data)

# 查看模型摘要
summary(model)

## 
## Call:
## lm(formula = LOGTPY ~ LOGNUMBED, data = filtered_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.87482 -0.02201  0.01517  0.05316  0.28862 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.17469    0.04537   -3.85  0.00014 ***
## LOGNUMBED    1.01923    0.01012  100.73  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09373 on 353 degrees of freedom
## Multiple R-squared:  0.9664, Adjusted R-squared:  0.9663 
## F-statistic: 1.015e+04 on 1 and 353 DF,  p-value: < 2.2e-16

# 提取回歸係數、t 統計量和 R-squared
coefficients <- summary(model)$coefficients
regression_coefficient <- coefficients["LOGNUMBED", "Estimate"]
t_statistic <- coefficients["LOGNUMBED", "t value"]
r_squared <- summary(model)$r.squared

線性回歸分析

我們建立了一個線性回歸模型，其中 LOGTPY 為因變數，LOGNUMBED 為自變數。回歸模型的主要結果如下：

• 回歸係數：
  • 截距項（Intercept）：-0.17469，該值表示當 LOGNUMBED 為 0 時，LOGTPY 的預測值。
  • LOGNUMBED 的回歸係數為 1.01923，這表明當 LOGNUMBED 增加 1 單位時，LOGTPY 預計會增加約 1.01923 單位。
• t 統計量：
  • 截距項的 t 值為 -3.85，對應的 p 值為 0.00014，顯示截距項在統計上顯著。
  • LOGNUMBED 的 t 值為 100.73，對應的 p 值小於 2e-16，這表明 LOGNUMBED 與 LOGTPY 之間存在非常顯著的正相關。
• 決定係數：
  • 決定係數（\(R^2\)）為 0.9664，表示大約 96.64% 的 LOGTPY 變異可以通過 LOGNUMBED 來解釋，這表明模型擬合效果非常好。

結果總結

從模型的結果來看，LOGNUMBED 對 LOGTPY 具有顯著的影響。回歸分析顯示，當 LOGNUMBED 增加時，LOGTPY 隨之增加。回歸係數為 1.01923，並且 t 統計量的結果顯示，這一關係在統計上極為顯著（t = 100.73，p-value < 2e-16）。此外，決定係數 \(R^2 = 0.9664\) 表明該模型能夠很好地解釋 LOGTPY 的變異。

這個回歸模型顯示了 LOGNUMBED 與LOGTPY之間強烈且顯著的正向關聯，並且模型擬合度非常高。

Test the following hypotheses at the 5 level of significance using a t -statistic. Also compute the corresponding p-value.

假設檢定

在這裡，我們將對回歸模型中的回歸係數 \(\beta_1\) 進行四個不同的假設檢定，這些檢定將幫助我們了解 \(\beta_1\) 是否顯著不同於某些值。

假設檢定 1： \(H_0 : \beta_1 = 0\) 對立假設 \(H_a : \beta_1 \neq 0\)

我們的零假設是 \(\beta_1 = 0\)，對立假設是 \(\beta_1 \neq 0\)。這是一個雙尾檢定，計算的 t 統計量如下：

\[ t = \frac{\hat{\beta_1} - 0}{SE(\hat{\beta_1})} \]

其中 \(\hat{\beta_1}\) 是估計值， \(SE(\hat{\beta_1})\) 是標準誤差。

# 假設回歸模型已經擬合並存儲在 'model' 中
# 提取回歸係數的估計值和標準誤差
beta1_hat <- summary(model)$coefficients[2, 1]
SE_beta1_hat <- summary(model)$coefficients[2, 2]

# (c i) 測試 H0: β1 = 0 對 Ha: β1 ≠ 0
t_stat_1 <- (beta1_hat - 0) / SE_beta1_hat
p_value_1 <- 2 * (1 - pt(abs(t_stat_1), df = df.residual(model)))
# 顯示結果
cat("c(i) p-value (H0: β1 = 0 vs Ha: β1 != 0): ", p_value_1, "\n")

## c(i) p-value (H0: β1 = 0 vs Ha: β1 != 0):  0

p_value_ci <- 0
result_ci <- ifelse(p_value_ci < 0.05, "Reject H0", "Fail to reject H0")
result_ci

## [1] "Reject H0"

假設檢定 2： \(H_0 : \beta_1 = 1\) 對立假設 \(H_a : \beta_1 \neq 1\)

在這個假設檢定中，我們的零假設是 \(\beta_1 = 1\)，對立假設是 \(\beta_1 \neq 1\)。這也是一個雙尾檢定，t 統計量公式類似：

\[ t = \frac{\hat{\beta_1} - 1}{SE(\hat{\beta_1})} \]

# (c ii) 測試 H0: β1 = 1 對 Ha: β1 ≠ 1
t_stat_2 <- (beta1_hat - 1) / SE_beta1_hat
p_value_2 <- 2 * (1 - pt(abs(t_stat_2), df = df.residual(model)))
cat("c(ii) p-value (H0: β1 = 1 vs Ha: β1 != 1): ", p_value_2, "\n")

## c(ii) p-value (H0: β1 = 1 vs Ha: β1 != 1):  0.05817311

p_value_cii <- 0.05817311
result_cii <- ifelse(p_value_cii < 0.05, "Reject H0", "Fail to reject H0")
result_cii

## [1] "Fail to reject H0"

假設檢定 3： \(H_0 : \beta_1 = 1\) 對立假設 \(H_a : \beta_1 > 1\)

這是單尾檢定，對立假設是 \(\beta_1 > 1\)。我們計算 t 統計量並查找單尾 p 值。

# (c iii) 測試 H0: β1 = 1 對 Ha: β1 > 1 (單尾檢定)
t_stat_3 <- (beta1_hat - 1) / SE_beta1_hat
p_value_3 <- 1 - pt(t_stat_3, df = df.residual(model))
cat("c(iii) p-value (H0: β1 = 1 vs Ha: β1 > 1): ", p_value_3, "\n")

## c(iii) p-value (H0: β1 = 1 vs Ha: β1 > 1):  0.02908655

p_value_ciii <- 0.02908655
result_ciii <- ifelse(p_value_ciii < 0.05, "Reject H0", "Fail to reject H0")
result_ciii

## [1] "Reject H0"

假設檢定 4： \(H_0 : \beta_1 = 1\) 對立假設 \(H_a : \beta_1 < 1\)

這也是單尾檢定，對立假設是 \(\beta_1 < 1\)。同樣，我們會計算 t 統計量並查找單尾 p 值。

# (c iv) 測試 H0: β1 = 1 對 Ha: β1 < 1 (單尾檢定)
t_stat_4 <- (beta1_hat - 1) / SE_beta1_hat
p_value_4 <- pt(t_stat_4, df = df.residual(model))
cat("c(iv) p-value (H0: β1 = 1 vs Ha: β1 < 1): ", p_value_4, "\n")

## c(iv) p-value (H0: β1 = 1 vs Ha: β1 < 1):  0.9709134

p_value_civ <- 0.9709134
result_civ <- ifelse(p_value_civ < 0.05, "Reject H0", "Fail to reject H0")
result_civ

## [1] "Fail to reject H0"

results <- data.frame(
  Hypothesis = c("H0: β1 = 0 vs Ha: β1 != 0", 
                 "H0: β1 = 1 vs Ha: β1 != 1", 
                 "H0: β1 = 1 vs Ha: β1 > 1", 
                 "H0: β1 = 1 vs Ha: β1 < 1"),
  p_value = c(0, 0.05817311, 0.02908655, 0.9709134),
  Decision = c(result_ci, result_cii, result_ciii, result_civ)
)
knitr::kable(results, caption = "Summary of Hypothesis Testing Results")

d(i). Suppose that there is a marginal change in LOGNUMBED of 2.

在回歸模型中，\(\hat{\beta_1} = 1.01923\) 是 \(\text{LOGNUMBED}\) 的回歸係數。這代表 \(\text{LOGNUMBED}\) 每變化一單位，\(\text{LOGTPY}\) 預期會改變 \(\hat{\beta_1}\) 單位。

當 \(\Delta \text{LOGNUMBED} = 2\) 時，對應的 \(\Delta \text{LOGTPY}\) 的點估計值可由以下公式計算：

\[ \Delta \text{LOGTPY} = \hat{\beta_1} \times \Delta \text{LOGNUMBED} \]

# 設定回歸係數
beta_1_hat <- 1.01923

# 設定 LOGNUMBED 的變化
delta_lognumbed <- 2

# 計算 LOGTPY 的變化
delta_logtpy <- beta_1_hat * delta_lognumbed

# 輸出結果
cat("當 ΔLOGNUMBED =", delta_lognumbed, "時，預期的 ΔLOGTPY 為：", delta_logtpy, "\n")

## 當 ΔLOGNUMBED = 2 時，預期的 ΔLOGTPY 為： 2.03846

d(ii). Provide a 95% confidence interval corresponding to the point estimate in part d(i).

# 計算標準誤
se_delta_logtpy <- summary(model)$coefficients["LOGNUMBED", "Std. Error"] 

# 計算 t 臨界值
alpha <- 0.05
t_critical_95 <- qt(1 - alpha / 2, df = df.residual(model))

# 計算 95% 置信區間
lower_bound_95 <- 1.01923 - t_critical_95 * se_delta_logtpy
upper_bound_95 <- 1.01923 + t_critical_95 * se_delta_logtpy

cat(lower_bound_95*2, upper_bound_95*2)

## 1.99866 2.07826

計算

1. 標準誤的計算：
   • 標準誤（Standard Error）反映了回歸係數 \(\hat{\beta_1}\) 的估計不確定性。
   • 使用 summary(model)$coefficients["LOGNUMBED", "Std. Error"] 提取 \(\text{LOGNUMBED}\) 的標準誤。
2. t 臨界值的計算：
   • \(t_{\alpha/2}\) 是計算信賴區間的關鍵值。
   • 95% 信心水準對應 \(\alpha = 0.05\)，我們用自由度 \(df = \text{模型的自由度}\) 計算： \[ t_{\alpha/2} = qt(1 - \alpha / 2, df) \]
3. 信賴區間公式： \[ \text{信賴區間} = \hat{\beta_1} \pm t_{\alpha/2} \cdot SE(\hat{\beta_1}) \]
   • 下界：\(\hat{\beta_1} - t_{\alpha/2} \cdot SE(\hat{\beta_1})\)
   • 上界：\(\hat{\beta_1} + t_{\alpha/2} \cdot SE(\hat{\beta_1})\)
4. 擴展到 \(\Delta \text{LOGNUMBED} = 2\)：
   • 原回歸係數 \(\hat{\beta_1}\) 對應 \(\Delta \text{LOGNUMBED} = 1\)。
   • 將上下界各乘以 2，得到 \(\Delta \text{LOGNUMBED} = 2\) 時的信賴區間。

d(iii). Provide a 99% confidence interval corresponding to the point estimate in part d(i).

# 計算標準誤
se_delta_logtpy <- summary(model)$coefficients["LOGNUMBED", "Std. Error"] 

# 計算 t 臨界值
alpha <- 0.01
t_critical_99 <- qt(1 - alpha / 2, df = df.residual(model))

# 計算 95% 置信區間
lower_bound_99 <- 1.01923 - t_critical_99 * se_delta_logtpy
upper_bound_99 <- 1.01923 + t_critical_99 * se_delta_logtpy

cat(lower_bound_99*2, upper_bound_99*2)

## 1.98605 2.09087

e(i). Find the predicted value of LOGTPY.

filtered_data <- data %>% 
  filter(CRYEAR == 2001) %>% 
  select(TPY, NUMBED, SQRFOOT)

filtered_data$LOGTPY <- log(filtered_data$TPY)
filtered_data$LOGNUMBED <- log(filtered_data$NUMBED)

# 擬合線性回歸模型
model2 <- lm(LOGTPY ~ LOGNUMBED, data = filtered_data)

# 查看模型摘要
summary(model2)

## 
## Call:
## lm(formula = LOGTPY ~ LOGNUMBED, data = filtered_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.87482 -0.02201  0.01517  0.05316  0.28862 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.17469    0.04537   -3.85  0.00014 ***
## LOGNUMBED    1.01923    0.01012  100.73  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09373 on 353 degrees of freedom
## Multiple R-squared:  0.9664, Adjusted R-squared:  0.9663 
## F-statistic: 1.015e+04 on 1 and 353 DF,  p-value: < 2.2e-16

# 已知參數
beta_0 <- summary(model2)$coefficients["(Intercept)", "Estimate"]
beta_1 <- summary(model2)$coefficients["LOGNUMBED", "Estimate"]
LOGNUMBED_star <- log(100)  # 假設值

# 預測值計算
LOGTPY_pred <- beta_0 + beta_1 * LOGNUMBED_star
LOGTPY_pred

## [1] 4.519037

e(ii). Obtain the standard error of the prediction.

library(psych)

filtered_data <- data %>% 
  filter(CRYEAR == 2001) %>% 
  select(TPY, NUMBED, SQRFOOT)

filtered_data$LOGTPY <- log(filtered_data$TPY)
filtered_data$LOGNUMBED <- log(filtered_data$NUMBED)

data2 <- filtered_data %>%
select(LOGNUMBED, LOGTPY)

describe(data2)

##           vars   n mean   sd median trimmed  mad  min  max range  skew kurtosis
## LOGNUMBED    1 355 4.46 0.49    4.5    4.45 0.49 2.89 6.12  3.23 -0.07     0.55
## LOGTPY       2 355 4.37 0.51    4.4    4.37 0.49 2.51 6.09  3.58 -0.14     0.56
##             se
## LOGNUMBED 0.03
## LOGTPY    0.03

summary_data <- describe(data2)
# 已知數據
model_summary <- summary(model)
sigma_res <- model_summary$sigma   # 回歸模型的標準誤
n <- 355              # 樣本數
X_star <- log(100)    # 假設的預測點 (X^* = log(100))
X_bar  <-  summary_data$mean[summary_data$vars == 1]  # 自變數 NUMBED 的均值
sd_NUMBED <- summary_data$sd[summary_data$vars == 1]    # NUMBED 的標準差

# 計算總變異量 (sum of squares)
sum_squares <- (n-1) * (sd_NUMBED^2)

# 計算預測標準誤
SE_pred <- sqrt(sigma_res^2 * (1 + 1/n + ((X_star - X_bar)^2 / sum_squares)))
SE_pred

## [1] 0.09387851

描述

在這個分析中，我們將計算回歸模型中的預測標準誤差（Standard Error of the Prediction）。預測標準誤差衡量了基於回歸模型進行單一預測時的誤差範圍。它是用來估計預測結果的可靠性。

預測標準誤差的計算公式如下：

\[ SE(\hat{Y}^*) = \sqrt{\hat{\sigma}^2 \left( 1 + \frac{1}{n} + \frac{(X^* - \bar{X})^2}{\sum_{i=1}^n (X_i - \bar{X})^2} \right)} \]

• 其中：
  • .\(\hat{\sigma}^2\) 是模型的殘差變異數。
  • .\(X^*\) 是預測時使用的自變數值。
  • .\(\bar{X}\) 是自變數 \(X\) 的均值。
  • .\(n\) 是樣本數。

e(iii). Obtain a 95% prediction interval for your prediction.

# 計算 t 臨界值
alpha <- 0.05
t_critical_95 <- qt(1 - alpha / 2, df = df.residual(model))

# 計算 95% 置信區間
lower_bound_95 <- LOGTPY_pred - t_critical_95 * SE_pred
upper_bound_95 <- LOGTPY_pred + t_critical_95 * SE_pred

cat(lower_bound_95, upper_bound_95)

## 4.334405 4.703668

描述

在這個分析中，我們將使用 t 分佈來計算 t 臨界值，並基於回歸模型計算 95% 置信區間。置信區間提供了一個範圍，預測的結果有 95% 的機會落在這個範圍內。

步驟：

1. 計算 t 臨界值：使用 qt() 函數來計算對應於給定顯著性水平 (\(\alpha = 0.05\)) 的 t 臨界值。這個值將基於自由度（\(df\)）來計算，通常自由度是樣本數減去模型參數的數量。
2. 計算 95% 置信區間：使用以下公式來計算 95% 置信區間：

\[ CI = \hat{Y}^* \pm t_{\alpha/2} \cdot SE(\hat{Y}^*) \]

• 其中：
  • \(\hat{Y}^*\) 是預測值。
  • \(t_{\alpha/2}\) 是來自 t 分佈的臨界值，對應於 \(\alpha/2\)。
  • \(SE(\hat{Y}^*)\) 是預測的標準誤差。

e(iv). Convert the point prediction in part e(i) and the prediction interval obtained in part e(iii) into total person years (through exponentiation).

lower_bound_95_total <- exp(lower_bound_95)
upper_bound_95_total <- exp(upper_bound_95)
LOGTPY_pred_total <- exp(LOGTPY_pred)


# 顯示結果
cat("預測的TPY: ", LOGTPY_pred_total, "\n")

## 預測的TPY:  91.74716

cat("95% 信賴區間: [", lower_bound_95_total, ", ", upper_bound_95_total, "]")

## 95% 信賴區間: [ 76.27956 ,  110.3512 ]

將預測值轉換為總TPY：將 LOGTPY_pred 進行指數運算來還原到原始尺度。

e(v). Obtain a prediction interval as in part e(iv), corresponding to a 90% level (in lieu of 95%).

# 假設我們已經有回歸模型 'model' 和預測值 'LOGTPY_pred' 以及預測標準誤差 'SE_pred'

# 顯著性水平
alpha <- 0.10  # 90% 置信區間對應的顯著性水平

# 計算 t 臨界值
t_critical_90 <- qt(1 - alpha / 2, df = df.residual(model))

# 計算 90% 置信區間
lower_bound_90 <- LOGTPY_pred - t_critical_90 * SE_pred
upper_bound_90 <- LOGTPY_pred + t_critical_90 * SE_pred

# 將預測的對數值轉換回原始的總人年 (指數運算)
lower_bound_90_total <- exp(lower_bound_90)
upper_bound_90_total <- exp(upper_bound_90)
LOGTPY_pred_total <- exp(LOGTPY_pred)

# 顯示結果
cat("預測的總人年: ", LOGTPY_pred_total, "\n")

## 預測的總人年:  91.74716

cat("90% 置信區間: [", lower_bound_90_total, ", ", upper_bound_90_total, "]")

## 90% 置信區間: [ 78.58759 ,  107.1103 ]

回歸模型的假設

1. 線性關係假設 (Linearity)
   假設自變量與因變量之間存在線性關係，即因變量是自變量的線性函數。這意味著回歸模型假設回歸線能夠充分描述自變量和因變量之間的關係。

   • 解釋：
     如果實際關係是非線性的，回歸模型的預測將會偏差，因為線性模型無法正確捕捉變數之間的真實關係。

   
   

2. 獨立性假設 (Independence)
   假設每個觀察值（數據點）之間是相互獨立的。換句話說，因變量和自變量的觀察值不應該受到其他觀察值的影響。

   • 解釋：
     如果觀察值之間存在依賴性（如時間序列數據中的自相關），那麼這會違反回歸模型的假設，並可能導致錯誤的結果或估計。

   
   

3. 同方差性假設 (Homoscedasticity)
   假設誤差項的變異數是恆定的，也就是說，誤差項的分佈不隨著自變量的變化而改變。這意味著不同水平的自變量所對應的誤差應該有相同的變異數。

   • 解釋：
     如果誤差項的變異數隨自變量的變化而改變（稱為異方差性），那麼回歸模型的估計將不再是最有效的，且可能導致信賴區間和檢定結果的偏誤。

   
   

4. 正態性假設 (Normality)
   假設誤差項遵循正態分佈。這對回歸分析的參數估計並不那麼重要，但對於進行假設檢定（例如t檢定和F檢定）是必要的。

   • 解釋：
     如果誤差項不符合正態分佈，則回歸模型的檢定結果可能會受到影響，尤其是在樣本量較小的情況下，這會使得統計檢定無法正確反映真實情況。

   
   

5. 無多重共線性假設 (No Multicollinearity)
   假設自變量之間不應該存在高相關性。當自變量高度相關時（即多重共線性），模型無法確定哪些自變量對因變量有重要影響，會導致不穩定的係數估計。

   • 解釋：
     如果自變量之間高度相關，會導致回歸係數的標準誤過大，使得模型難以正確估計每個自變量對因變量的影響，從而降低模型的解釋能力。

   
   

6. 無遺漏變數假設 (No Omitted Variables)
   假設所有相關的自變量都被納入模型中。如果忽略了對因變量有重要影響的自變量，這會導致模型偏誤，並且估計結果不準確。

   • 解釋：
     遺漏變數會使得回歸模型的估計產生偏誤，因為它可能會將一些重要的變異來源錯誤地歸因於其他變量。

   
   

結論

這些假設是回歸分析的基礎，違反這些假設可能會導致模型的偏誤、效率低下或不準確的結論。在進行回歸分析時，我們需要檢查數據是否滿足這些假設，並根據需要進行數據處理或選擇更合適的分析方法。

i. 描述所選兩匹馬之間的因變量關係：

在每場比賽中選擇兩匹馬，這兩匹馬的因變量（例如比賽成績、名次）可能是相關的，而非相互獨立的。這種相關性可能來自於比賽條件（如天氣、賽道質量）或馬匹之間的競爭關係。例如，如果一匹馬表現良好，它可能會影響另一匹馬的表現，因為它們可能在相同的賽道上競爭，並受到相似的外部因素影響。

ii. 說明這如何違反回歸模型假設：

回歸模型的一個基本假設是數據點之間的獨立性，即每個觀察值（在這裡是每匹馬的表現）應該相互獨立。然而，當從同一場比賽中選擇兩匹馬時，這兩匹馬的表現很可能是相關的，因此違反了獨立性假設。這會影響回歸模型的結果，可能導致偏差的估計和錯誤的結論。

a. Compute descriptive statistics for TPY, NUMBED, and SQRFOOT.

# 使用 read_csv() 讀取 CSV 檔案
data <- read_csv("C:/Users/霍致光/Desktop/大學/應用機率模型/應用期末報告/HKHorse.csv")

# 使用 summary() 一次性分析
summary(data)

##      FINISH             WIN          
##  Min.   :0.00000   Min.   :0.002691  
##  1st Qu.:0.00000   1st Qu.:0.036856  
##  Median :0.00000   Median :0.079973  
##  Mean   :0.09514   Mean   :0.099565  
##  3rd Qu.:0.00000   3rd Qu.:0.133053  
##  Max.   :1.00000   Max.   :0.649985

# 使用 describe() 進行更詳細的描述性統計分析
describe(data)

##        vars   n mean   sd median trimmed  mad min  max range skew kurtosis   se
## FINISH    1 925  0.1 0.29   0.00    0.00 0.00   0 1.00  1.00 2.76      5.6 0.01
## WIN       2 925  0.1 0.08   0.08    0.09 0.07   0 0.65  0.65 1.63      4.0 0.00

Note that the standard deviation of FINISH is greater than that of WIN, even though the sample means are about the same. For the variable FINISH, what is the relationship between the sample mean and standard deviation?

1. FINISH 變數的特性

FINISH 是一個二元變數，只有兩個可能的值：0 和 1，分別代表馬匹沒有贏得比賽和贏得比賽。因此，FINISH 的樣本均值（\(\mu_y\)）表示的是馬匹贏得比賽的概率。

樣本均值的計算公式如下：

\[ \mu_y = \frac{\sum y_i}{n} \]

其中 \(y_i\) 為每個觀察值，\(n\) 為觀察值的數量。

2. 標準差與均值的關係

對於二元變數（0 或 1），標準差可以使用以下公式計算：

\[ \text{Standard Deviation} = \sqrt{\mu_y (1 - \mu_y)} \]

其中 \(\mu_y\) 是 FINISH 的樣本均值，即馬匹獲勝的概率。

3. 計算標準差

根據數據，FINISH 的樣本均值是 0.1。代入公式計算標準差：

\[ \text{Standard Deviation} = \sqrt{0.1 \times (1 - 0.1)} = \sqrt{0.1 \times 0.9} = \sqrt{0.09} = 0.3 \]

因此，FINISH 的標準差為 0.3。

4. FINISH 標準差大於 WIN 的標準差

儘管 FINISH 和 WIN 的樣本均值相似（約 0.1），但 FINISH 的標準差較大，這是由於 FINISH 變數的二元性質所致。FINISH 的值只能是 0 或 1，且這兩個值之間的差距較大，因此標準差較大。

而 WIN 變數是連續的，其取值範圍通常比較狹窄（例如在 0 到 0.65 之間），因此標準差較小。

5. 總結

• FINISH 是二元變數，標準差與均值之間有直接的數學關係，且當均值較小時，標準差通常較大。
• 儘管 FINISH 和 WIN 的樣本均值相似，但 FINISH 變數的標準差較大，這是由於其二元性質造成的。

Note that the sample mean is larger for horses that won (FINISH = 1) than for those that lost (FINISH = 0). Interpret this result.

計算 WIN 變數按 FINISH 分層的統計分析，並注意樣本均值對於贏馬（FINISH = 1）的樣本較大。解釋這一結果。

# 加載必要的包
library(tidyverse)

# 讀取數據
data <- read_csv("C:/Users/霍致光/Desktop/大學/應用機率模型/應用期末報告/HKHorse.csv")

# 按 FINISH 分組計算 WIN 的摘要統計
win_summary <- data %>%
  group_by(FINISH) %>%
  summarise(
    mean_win = mean(WIN, na.rm = TRUE),
    sd_win = sd(WIN, na.rm = TRUE),
    min_win = min(WIN, na.rm = TRUE),
    max_win = max(WIN, na.rm = TRUE),
    count = n()
  )

# 查看統計結果
print(win_summary)

## # A tibble: 2 × 6
##   FINISH mean_win sd_win min_win max_win count
##    <dbl>    <dbl>  <dbl>   <dbl>   <dbl> <int>
## 1      0   0.0931 0.0764 0.00269   0.469   837
## 2      1   0.161  0.104  0.00613   0.650    88

1. 均值解釋

• 當 FINISH = 1（贏馬）時，WIN 的平均值（0.161）高於當 FINISH = 0（未贏馬）時的平均值（0.0931）。
  
• 這表明人群對於贏馬馬匹的評估平均贏馬概率更高，顯示出人群預測與實際結果之間存在一定的正相關性。

2. 標準差解釋

• 贏馬的標準差（0.104）大於未贏馬的標準差（0.0764）。
  
• 這說明對於贏馬的馬匹，人群的預測分佈更分散，可能因為部分馬匹的表現超出了預期。例如，低評估的馬匹贏馬，或高評估的馬匹未能贏馬。

3. 極端值觀察

• 當 FINISH = 1 時，WIN 的最小值為 0.00613，表示某些贏馬的馬匹事前被嚴重低估。
  
• 當 FINISH = 0 時，WIN 的最大值為 0.469，表明即使是未贏馬的馬匹，人群中也有一些對其贏馬的信心較高。

4. 總結與結論

• 人群對馬匹的贏馬概率預測具備一定準確性，因為贏馬馬匹的平均預測值高於未贏馬馬匹。
  
• 然而，預測的準確性存在局限，部分結果顯示人群對贏馬概率的評估受隨機性影響，如比賽條件和馬匹狀態的變化。

i. Is WIN a statistically significant predictor of FINISH?

使用 WIN 預測 FINISH，估計線性概率模型，並檢查 WIN 是否為 FINISH 的統計顯著預測變量。

# 加載必要的包
library(tidyverse)
library(broom)

# 讀取數據
data <- read_csv("C:/Users/霍致光/Desktop/大學/應用機率模型/應用期末報告/HKHorse.csv")

# 建立線性概率模型
lpm_model <- lm(FINISH ~ WIN, data = data)

summary(lpm_model)

## 
## Call:
## lm(formula = FINISH ~ WIN, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41878 -0.11175 -0.06672 -0.03128  0.98661 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.008022   0.014747   0.544    0.587    
## WIN         0.874936   0.114408   7.647 5.14e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2848 on 923 degrees of freedom
## Multiple R-squared:  0.05959,    Adjusted R-squared:  0.05857 
## F-statistic: 58.48 on 1 and 923 DF,  p-value: 5.14e-14

# 提取回歸係數及顯著性檢定
tidy(lpm_model)

## # A tibble: 2 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)  0.00802    0.0147     0.544 5.87e- 1
## 2 WIN          0.875      0.114      7.65  5.14e-14

ii. How well does this model fit the data using the usual goodnessof-fit statistic?

# 計算模型的 R²
r_squared <- summary(lpm_model)$r.squared

# 輸出結果
cat("模型的 R² 值為：", r_squared)

## 模型的 R² 值為： 0.05958747

iii. For this estimated model, is it possible for the fitted values to lie outside the interval [0, 1]?

e. Estimate a logistic regression model, using WIN to predict FINISH.Is WIN a statistically significant predictor of FINISH?

使用 WIN 作為自變量，FINISH 作為因變量，建立邏輯回歸模型：

# 加載所需的套件
library(broom)

# 建立邏輯回歸模型
logit_model <- glm(FINISH ~ WIN, data = data, family = binomial)

# 查看模型摘要
summary(logit_model)

## 
## Call:
## glm(formula = FINISH ~ WIN, family = binomial, data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.1840     0.2004 -15.885  < 2e-16 ***
## WIN           7.6467     1.1470   6.667 2.61e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 581.38  on 924  degrees of freedom
## Residual deviance: 537.50  on 923  degrees of freedom
## AIC: 541.5
## 
## Number of Fisher Scoring iterations: 5

# 提取模型結果
tidy_logit <- tidy(logit_model)

# 顯示 WIN 的統計結果
tidy_logit

## # A tibble: 2 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    -3.18     0.200    -15.9  8.11e-57
## 2 WIN             7.65     1.15       6.67 2.61e-11

i.For each model, provide fitted values at WIN = 0, 0.01, 0.05, 0.10,and 1.0.

# 線性概率模型的擬合值計算
linear_fitted <- function(WIN) {
  intercept <- 0.00802
  coefficient <- 0.875
  return(intercept + coefficient * WIN)
}

# 計算擬合值
WIN_values <- c(0, 0.01, 0.05, 0.10, 1.0)
linear_fitted_values <- sapply(WIN_values, linear_fitted)
linear_fitted_values

## [1] 0.00802 0.01677 0.05177 0.09552 0.88302

# 邏輯回歸模型的擬合值計算

logit_fitted <- function(WIN) {
  intercept <- -3.18
  coefficient <- 7.65
  return(1 / (1 + exp(-(intercept + coefficient * WIN))))
}

# 計算擬合值
logit_fitted_values <- sapply(WIN_values, logit_fitted)
logit_fitted_values

## [1] 0.03992533 0.04296311 0.05745942 0.08203600 0.98868224

ii. Plot fitted values from the linear probability model versus fitted values from the logistic regression model.

# 繪製擬合值比較圖
plot(WIN_values, linear_fitted_values, col = "blue", pch = 16, 
     xlab = "WIN", ylab = "Fitted Values", 
     main = "Comparison of Fitted Values: LPM vs Logit")
points(WIN_values, logit_fitted_values, col = "red", pch = 17)
legend("bottomright", legend = c("Linear Probability Model", "Logistic Regression Model"), 
       col = c("blue", "red"), pch = c(16, 17))

i. Plot the difference FINISH – WIN versus WIN.

# 計算 FINISH 和 WIN 之間的差異
data$Difference <- data$FINISH - data$WIN

# 繪製差異與 WIN 的關係圖
library(ggplot2)

ggplot(data, aes(x = WIN, y = Difference)) +
  geom_point(alpha = 0.5) +
  labs(title = "Difference between FINISH and WIN vs. WIN",
       x = "Crowd's Prior Probability (WIN)",
       y = "Difference (FINISH - WIN)") +
  theme_minimal()

1. 對於差異較大的馬匹進行投注

• 如果 FINISH - WIN 的差異接近 1，這表示群眾對該匹馬的勝率預測較低，而實際上該匹馬可能有較高的勝算。在這種情況下，您可以選擇投注這些馬匹，因為這樣的馬匹往往會有較高的賠率，從而提供較高的潛在回報。
• 策略建議：選擇 FINISH - WIN 值大於某一閾值的馬匹進行投注。例如，當 FINISH - WIN > 0.5 時，可以考慮投注這些馬匹。

2. 對於差異較小的馬匹進行投注

• 如果 FINISH - WIN 的差異接近 0，這表示群眾的預測與實際結果相符，馬匹的勝算較為穩定。這些馬匹的賠率通常較低，但投注這些馬匹的風險較小，適合尋求穩定回報的投注者。
• 策略建議：選擇 FINISH - WIN 值接近 0 的馬匹進行投注，特別是在賠率相對較低的情況下，這樣可以減少風險並獲得穩定的回報。

3. 根據群眾預測進行反向投注

• 當 FINISH - WIN 差異較大時，您也可以考慮反向投注，即在群眾預測勝利機率較高時，選擇投注那些實際結果不太可能勝出的馬匹。這樣做可能會有更高的賠率，但同時風險較高。
• 策略建議：選擇 WIN 值較高的馬匹進行反向投注，尤其是在 FINISH - WIN 差異較大的情況下。這樣的馬匹可能會讓您以較高的賠率賺取潛在回報。

4. 風險控制

• 在採取上述策略時，風險控制是非常重要的。投注者應該根據自己的風險承受能力來決定投注金額，避免過度集中於某一種策略。
• 策略建議：根據 FINISH - WIN 的分佈情況來分散投注，避免將所有資金投入於差異較大的馬匹，這樣可以減少單一結果的風險。

結論

• 根據 FINISH - WIN 差異，可以設計出多種策略來進行賽馬投注。最關鍵的在於了解這些差異背後的含義，並且能夠靈活調整策略，以最大化回報並控制風險。投注者應根據自己的風險偏好來選擇合適的策略，並且不斷根據賽事結果進行調整。