```{r}
# ЧАСТЬ 5. Оценка порядковой пробит-модели
model <- polr(as.factor(sub_type) ~ income + age + internet + series + health + male +
              marriage + residence + cat + news + TV + time, data = data, Hess = TRUE)

# Посмотрим на результаты
print(model)
```
```{r}
# Выбор случайных значений для каждого признака
individual <- data[sample(nrow(data), 1), -c(11,12)]

# Прогнозирование вероятностей
predicted_probs <- predict(model, newdata=individual, type="probs")

# Вероятность для обычной подписки
prob_ordinary <- predicted_probs[2]
print(prob_ordinary)
```


```{r}
# Оценка порядковой пробит-модели
model <- polr(as.factor(sub_type) ~ income + age + internet + series + health + male +
    marriage + residence + cat + news + TV + time, data=data, Hess=TRUE)
```


```{r}
# Коэффициенты из вашей порядковой пробит-модели
coef <- c(income = 2.084015e-05, age = 1.471304e-02, internet = 1.627123e+00, series = 2.306277e-01,
          healthgood = 5.127466e-01, healthmedium = 3.337718e-02, male = 3.443925e+00, marriage = -4.772577e-02,
          residenceCity = -6.251955e-01, residenceVillage = -2.100985e+00, cat = -9.918646e-03,
          news = 7.752378e-02, TV = -2.171019e+00, time = -4.459438e-04)

# Предельные эффекты
effects <- coef / (1 + exp(-coef))

# Вывод результатов
effects
```

```{r}
# Создание дамми-переменных
data$healthgood <- as.numeric(data$health == "good")
data$healthmedium <- as.numeric(data$health == "medium")
data$residenceCity <- as.numeric(data$residence == "City")
data$residenceVillage <- as.numeric(data$residence == "Village")

# Выбор случайного индивида из данных
random_individual <- data[sample(nrow(data), 1),]

# Коэффициенты из вашей порядковой пробит-модели
coef <- c(income = 2.084015e-05, age = 1.471304e-02, internet = 1.627123e+00, series = 2.306277e-01,
          healthgood = 5.127466e-01, healthmedium = 3.337718e-02, male = 3.443925e+00, marriage = -4.772577e-02,
          residenceCity = -6.251955e-01, residenceVillage = -2.100985e+00, cat = -9.918646e-03,
          news = 7.752378e-02, TV = -2.171019e+00, time = -4.459438e-04)

# Предельные эффекты для непрерывных переменных
continuous_effects <- coef / (1 + exp(-coef))

# Предельные эффекты для дискретных переменных
discrete_vars <- c("male", "marriage", "residenceCity", "residenceVillage", "cat", "TV", "healthgood", "healthmedium")
discrete_effects <- numeric(length(discrete_vars))
for (i in 1:length(discrete_vars)) {
    x0 <- as.numeric(random_individual[, names(coef)])
    x1 <- x0
    x1[which(names(coef) == discrete_vars[i])] <- 1 - x0[which(names(coef) == discrete_vars[i])]
    prob_x0 <- pnorm(sum(coef * x0))
    prob_x1 <- pnorm(sum(coef * x1))
    discrete_effects[i] <- prob_x1 - prob_x0
}
names(discrete_effects) <- discrete_vars

# Вывод результатов
list(Continuous = continuous_effects, Discrete = discrete_effects)
```
```{r}
# Полная модель
full_model <- polr(as.factor(sub_type) ~ income + age + internet + series + health + male + marriage + residence + cat + news + TV + time, data = data)

# Ограниченная модель (предположим, мы хотим сравнить коэффициенты для переменных 'income' и 'age')
data$combined_var <- data$income + data$age
restricted_model <- polr(as.factor(sub_type) ~ combined_var + internet + series + health + male + marriage + residence + cat + news + TV + time, data = data)

# LR статистика
lr_stat <- -2 * (logLik(restricted_model) - logLik(full_model))

# P-значение
p_value <- 1 - pchisq(lr_stat, df = 1)

print(lr_stat)
```
```{r}
print(p_value)
```
<!-- ЧАСТЬ 6 -->
```{r}
library(nnet)
# Оценка мультиномиальной логит-модели
mnl_model <- multinom(sub_type ~ income + age + internet + series + health + male + marriage + residence + cat + news + TV + time, data = data)

# Вывод результатов
print(mnl_model)
```
```{r}

# Выбор случайных значений для каждого признака
individual <- data[sample(nrow(data), 1),]
print(individual)
# Применение модели к данным индивида
predicted_probs <- predict(mnl_model, newdata = individual, type = "probs")

# Вероятность выбора обычной подписки
ans <- predicted_probs[2]

print(ans)

```

```{r}
# Выбор случайного индивида из данных
random_individual <- data[sample(nrow(data), 1),]

# Предельные эффекты для непрерывных переменных
continuous_vars <- c("income", "age", "internet", "series")
continuous_effects <- numeric(length(continuous_vars))
predicted_probs <- predict(mnl_model, newdata = random_individual, type = "probs")
print(names(coef))
for (i in 1:length(continuous_vars)) {
    beta_k <- coef[continuous_vars[i]] # обновленное извлечение коэффициентов
    prob_k <- predicted_probs[2]
    prob_0 <- predicted_probs[1]
    continuous_effects[i] <- (beta_k * prob_k) - (0 * prob_0) # 0 для базовой категории
}

# Предельные эффекты для дискретных переменных
discrete_vars <- c("male", "marriage", "residenceCity", "residenceVillage", "cat", "TV", "healthgood", "healthmedium")
discrete_effects <- numeric(length(discrete_vars))
for (i in 1:length(discrete_vars)) {
    x0 <- as.numeric(random_individual[, discrete_vars[i]])
    x1 <- 1 - x0
    random_individual[, discrete_vars[i]] <- x1
    prob_x1 <- predict(mnl_model, newdata = random_individual, type = "probs")[2]
    discrete_effects[i] <- prob_x1 - prob_k
}

# Вывод результатов
list(Continuous = setNames(continuous_effects, continuous_vars),
     Discrete = setNames(discrete_effects, discrete_vars))
```


```{r}
install.packages("sampleSelection")

# Оценка модели Хекмана
result <- heckit(subscription ~ age + income + education,  # Уравнение участия
                 watch_time ~ age + income,                # Уравнение результата
                 data = your_data_frame)

# Печать результатов
summary(result)
```