"Interation Model" = log_interaction_lm
),
stars = TRUE,
fmt = 5,
title = 'Comparsion Table between Two Models',
output = "html")
# we add a new column in sales dataframe
sales$predicted <- predict(log_interaction_lm, data = sales)
# check the prediction
sales$predicted[1]
sales$predicted_exp <- exp(sales$predicted)
sales$predicted_exp[1]
sales$condition
sales$condition %>% unique()
AIC(simple_lm)
cor(sales[, c("sale_price", "sqft", "year_built", "condition")])
sales_long <- sales |>
pivot_longer(
cols = c(sqft, year_built),
names_to = "variable",
values_to = "xvalue"
)
ggplot(sales_long, aes(x = xvalue, y = sale_price)) +
geom_point(size = 0.3, alpha = 0.5, color = "#e8e3d3") +
geom_smooth(linewidth = 0.5, method = "lm", formula = "y~x", color = "#4b2e83") +
facet_wrap(~ variable, scales = "free_x") +
labs(
x = "Explanatory Variable",
y = "Sale Price ($)",
title = "Sale Price vs. Square Footage of House (Sqft) and Year Built"
) +
theme_minimal()
sales_long <- sales |>
pivot_longer(
cols = c(sqft, year_built, condition),
names_to = "variable",
values_to = "xvalue"
)
ggplot(sales_long, aes(x = xvalue, y = sale_price)) +
geom_point(size = 0.3, alpha = 0.5, color = "#e8e3d3") +
geom_smooth(linewidth = 0.5, method = "lm", formula = "y~x", color = "#4b2e83") +
facet_wrap(~ variable, scales = "free_x") +
labs(
x = "Explanatory Variable",
y = "Sale Price ($)",
title = "Sale Price vs. Square Footage of House (Sqft) and Year Built"
) +
theme_minimal()
sales_long <- sales |>
pivot_longer(
cols = c(sqft, year_built),
names_to = "variable",
values_to = "xvalue"
)
ggplot(sales_long, aes(x = xvalue, y = sale_price)) +
geom_point(size = 0.3, alpha = 0.5, color = "#e8e3d3") +
geom_smooth(linewidth = 0.5, method = "lm", formula = "y~x", color = "#4b2e83") +
facet_wrap(~ variable, scales = "free_x") +
labs(
x = "Explanatory Variable",
y = "Sale Price ($)",
title = "Sale Price vs. Square Footage of House (Sqft) and Year Built"
) +
theme_minimal()
# lm means linear model
# we don't need to add intercept, R will add it for us
simple_lm <- lm(data = sales, sale_price ~ sqft + year_built)
summary(simple_lm)
confint(simple_lm)
example_x <- seq(1, 100, by = 1)
example_y <- 3 + 2 * example_x + rnorm(100, mean = 0, sd = 10)
model <- lm(example_y ~ example_x)
summary(model)
plot(model)
example_x <- seq(1, 100, by = 1)
example_y <- 3 + 2 * example_x + rnorm(100, mean = 0, sd = 2)
example_model <- lm(example_y ~ example_x)
summary(model)
plot(model)
example_x <- seq(1, 100, by = 1)
example_y <- 3 + 2 * example_x + rnorm(100, mean = 0, sd = 2)
example_model <- lm(example_y ~ example_x)
plot(example_model)
example_x <- seq(1, 100, by = 1)
example_y <- 3 + 2 * example_x + rnorm(100, mean = 0, sd = 2)
example_model <- lm(example_y ~ example_x)
plot(example_model, which = 1)
plot(simple_lm, which = 1)
plot(example_model, which = 1, title("Perfect"))
plot(example_model, which = 1, title("Perfect Model"))
plot(simple_lm, which = 1, title("Our Linear Model"))
plot(example_model, which = 2, title("Perfect Model"))
plot(simple_lm, which = 2, title("Our Linear Model"))
plot(example_model, which = 5, title("Perfect Model"))
plot(simple_lm, which = 5, title("Our Linear Model"))
#install.packages('car')
library(car)
vif(simple_lm)
coef(log_simple_lm)
summary(log_simple_lm)
log_simple_lm <- lm(data = sales, log(sale_price) ~ sqft + year_built)
plot(log_simple_lm)
knitr::opts_chunk$set(echo = TRUE)
hist(sales$sale_price)
hist(log(sales$sale_price))
dat <- data.frame(responses = c(2600, 2900, 2000, 2200,3200, 35000,23000,20000,
30000,27000,2900000,2300000,1700000,2900000,2000000),
treatment = as.factor(c(rep("1",5),rep("2",5),rep("3",5))))
library(ggplot2)
View(dat)
ggplot(data = dat) +
geo_point()
ggplot(data = dat) +
geom_point()
ggplot(data = dat, aes(treatment, responses)) +
geom_point()
ggplot(data = dat, aes(treatment, responses)) +
geom_point()
ggplot(data = dat, aes(treatment, responses)) +
geom_point() +
facet_grid()
ggplot(data = dat, aes(treatment, responses)) +
geom_point() +
facet_grid(treatment)
ggplot(data = dat, aes(treatment, responses)) +
geom_point() +
facet_grid(. ~ treatment)
ggplot(data = dat, aes(treatment, responses)) +
geom_point() +
facet_wrap(~ treatment)
ggplot(dat, aes(x = treatment, y = responses)) +
geom_point() +
facet_wrap(~ treatment, scales = "free_y") +
theme_bw()
fit <- aov(responses ~ treatment, data = dat)
plot(fit, which = 2)
plot(fit, which = 1)
fit <- aov(responses ~ treatment, data = dat)
plot(fit)
summary(fit)
fit <- aov(responses ~ treatment, data = dat)
summary(fit)
plot(fit)
library(MASS)
bc <- boxcox(fit)
lambda <- bc$x[which.max(bc$y)]
lambda
library(MASS)
bc <- boxcox(fit)
lambda <- bc$x[which.max(bc$y)]
lambda
means <- tapply(y, group, mean)
means <- tapply(dat$responses, dat$treatment, mean)
sds   <- tapply(dat$responses, dat$treatment, sd)
plot(log(sds) ~ log(means), pch=19)
means <- tapply(dat$responses, dat$treatment, mean)
sds   <- tapply(dat$responses, dat$treatment, sd)
plot(log(sds) ~ log(means), pch=19)
fit_alpha <- lm(log(sds) ~ log(means))
summary(fit_alpha)
means <- tapply(dat$responses, dat$treatment, mean)
sds   <- tapply(dat$responses, dat$treatment, sd)
fit_alpha <- lm(log(sds) ~ log(means))
summary(fit_alpha)
means <- tapply(dat$responses, dat$treatment, mean)
sds   <- tapply(dat$responses, dat$treatment, sd)
fit_alpha <- lm(log(sds) ~ log(means))
summary(fit_alpha)
dat$responses_log <- log(dat$responses)
dat$responses_log <- log(dat$responses)
fit2 <- aov(responses_log ~ treatment, data = dat)
summary(fit2)
plot(fit2)
load('lentil.dat')
d.len <- read.table("lentil.dat", header = TRUE)
dat <- read.table("lentil.dat", header = TRUE)
View(dat)
d.len <- read.table("lentil.dat", header = TRUE)
d.len$TR <- factor(d.len$TR)
View(d.len)
fit <- aov(Y ~ TR, data = d.len)
summary(fit.len)
summary(fit)
plot(fit)
plot(fit, which = c(1,2))
hist(residuals(fit.len))
hist(residuals(fit))
hist(residuals(fit))
plot(fit, which = c(1,2))
# Define contrast matrix
mat.contr <- rbind(
C1 = c(-6, +1, +1, +1, +1, +1, +1),
C2 = c( 0, -1, -1, -1, +1, +1, +1),
C3 = c( 0, +2, -1, -1, +2, -1, -1),
C4 = c( 0,  0, -1, +1,  0, -1, +1),
C5 = c( 0, -2, +1, +1, +1, -2, +1),
C6 = c( 0,  0, +1, -1,  0, -1, +1)
)
# Check orthogonality
t(mat.contr) %*% mat.contr
# Define contrast matrix
mat.contr <- rbind(
C1 = c(-6, +1, +1, +1, +1, +1, +1),
C2 = c( 0, -1, -1, -1, +1, +1, +1),
C3 = c( 0, +2, -1, -1, +2, -1, -1),
C4 = c( 0,  0, -1, +1,  0, -1, +1),
C5 = c( 0, -2, +1, +1, +1, -2, +1),
C6 = c( 0,  0, +1, -1,  0, -1, +1)
)
# Check orthogonality
t(mat.contr) %*% mat.contr
# Define contrast matrix
mat <- rbind(
C1 = c(-6, +1, +1, +1, +1, +1, +1),
C2 = c( 0, -1, -1, -1, +1, +1, +1),
C3 = c( 0, +2, -1, -1, +2, -1, -1),
C4 = c( 0,  0, -1, +1,  0, -1, +1),
C5 = c( 0, -2, +1, +1, +2, -1, -1),
C6 = c( 0,  0, +1, -1,  0, -1, +1)
)
# Check orthogonality
t(mat.contr) %*% mat.contr
# Check orthogonality
t(mat) %*% mat
# Check orthogonality
mat.contr %*% t(mat.contr)
# Define contrast matrix
mat <- rbind(
C1 = c(-6, +1, +1, +1, +1, +1, +1),
C2 = c( 0, -1, -1, -1, +1, +1, +1),
C3 = c( 0, +2, -1, -1, +2, -1, -1),
C4 = c( 0,  0, -1, +1,  0, -1, +1),
C5 = c( 0, -2, +1, +1, +2, -1, -1),
C6 = c( 0,  0, +1, -1,  0, -1, +1)
)
# Check orthogonality
mat.contr %*% t(mat.contr)
mat.contr <- rbind(
C1 = c(-6, +1, +1, +1, +1, +1, +1),
C2 = c( 0, -1, -1, -1, +1, +1, +1),
C3 = c( 0, +2, -1, -1, +2, -1, -1),
C4 = c( 0,  0, -1, +1,  0, -1, +1),
C5 = c( 0, -2, +1, +1, +2, -1, -1),
C6 = c( 0,  0, +1, -1,  0, -1, +1)
)
mat.contr %*% t(mat.contr)
mat_contrasts <- rbind(
C1 = c(-6, +1, +1, +1, +1, +1, +1),
C2 = c( 0, -1, -1, -1, +1, +1, +1),
C3 = c( 0, +2, -1, -1, +2, -1, -1),
C4 = c( 0,  0, -1, +1,  0, -1, +1),
C5 = c( 0, -2, +1, +1, +2, -1, -1),
C6 = c( 0,  0, +1, -1,  0, -1, +1)
)
mat_contrasts %*% t(mat_contrasts)
fit.mc <- glht(fit.len, linfct = mcp(TR = mat.contr))
library(multcomp)
install.packages("multcomp")
library(multcomp)
fit.mc <- glht(fit.len, linfct = mcp(TR = mat.contr))
library(multcomp)
fit.mc <- glht(fit.len, linfct = mcp(TR = mat_contrasts))
fit.mc <- glht(fit, linfct = mcp(TR = mat_contrasts))
library(multcomp)
fit.mc <- glht(fit, linfct = mcp(TR = mat_contrasts))
summary(fit.mc, test = adjusted("none"))
summary(fit.mc, test = adjusted("bonferroni"))
summary(fit.mc, test = adjusted("bonferroni"))
library(multcomp)
fit.mc <- glht(fit, linfct = mcp(TR = mat_contrasts))
summary(fit.mc, test = adjusted("bonferroni"))
y <- c(9, 12, 10, 8, 15, 20, 21, 23, 17, 30, 6, 5, 8, 16, 7)
type <- c(...)
y <- c(9, 12, 10, 8, 15, 20, 21, 23, 17, 30, 6, 5, 8, 16, 7)
type <-  c(rep(1,5), rep(2,5), rep(3,5))
circ <- data.frame(Type = type, Y = y)
circ$Type <- as.factor(circ$Type)
fit <- aov(Y ~ Type, data = circ)
summary(fit)
plot(fit)
fit <- aov(Y ~ Type, data = circ)
summary(fit)
par(mfrow=c(2,2))
plot(fit)
summary(fit)
fit <- aov(Y ~ Type, data = circ)
summary(fit)
plot(fit)
plot(fit, which = c(1,2))
mat_contrasts <- rbind(
C1 = c(-1, +2, -1),
C2 = c(+1, +0, -1)
)
mat_contrasts %*% t(mat_contrasts)
circ.mc <- glht(fit, linfct = mcp(Type = mat_contrasts))
summary(circ.mc, test = adjusted("holm"))
mat_contrasts <- rbind(
C1 = c(-1, +2, -1), # this compare type 2 with the remaining
C2 = c(+1, +0, -1),
C3 = c(+0, +1, -1)
)
mat_contrasts %*% t(mat_contrasts)
mat_contrasts <- rbind(
C1 = c(-1, +2, -1), # this compare type 2 with the remaining
C2 = c(+1, +0, -1))
mat_contrasts %*% t(mat_contrasts)
circ.mc <- glht(fit, linfct = mcp(Type = mat_contrasts))
summary(circ.mc, test = adjusted("holm"))
#install.packages("tidymodels") # you only need to install once
library(tidymodels)
# 8. 汇总表格（平均性能） -------------------------------------------
collect_metrics(results) %>%
select(wflow_id, mean, std_err, .metric) %>%
pivot_wider(names_from = .metric, values_from = mean)
sales_pred <- sales %>%
mutate(view_mountains = view_rainier + view_olympics + view_cascades + view_territorial,
view_water = view_sound + view_lakewash + view_lakesamm + view_otherwater)
# remove non-predictive or unique identifiers
sales_pred <- sales_pred %>%
select(-sale_id, -pinx, -sale_date, -sale_warning, -zoning, -join_status, -join_year, -longitude, -latitude, -city, -area, -land_val, -imp_val, -sqft_1, -sqft_fbsmt, -submarket, -view_sound, -view_lakewash, -view_lakesamm, -view_otherwater, -view_rainier, -view_olympics, -view_cascades, -view_territorial, -view_other, -subdivision, -predicted, -predicted_exp)
sales_pred <- sales %>%
mutate(view_mountains = view_rainier + view_olympics + view_cascades + view_territorial,
view_water = view_sound + view_lakewash + view_lakesamm + view_otherwater)
# remove non-predictive or unique identifiers
sales_pred <- sales_pred %>%
select(-sale_id, -pinx, -sale_date, -sale_warning, -zoning, -join_status, -join_year, -longitude, -latitude, -city, -area, -land_val, -imp_val, -sqft_1, -sqft_fbsmt, -submarket, -view_sound, -view_lakewash, -view_lakesamm, -view_otherwater, -view_rainier, -view_olympics, -view_cascades, -view_territorial, -view_other, -subdivision, -predicted, -predicted_exp)
sales_pred <- sales %>%
mutate(view_mountains = view_rainier + view_olympics + view_cascades + view_territorial,
view_water = view_sound + view_lakewash + view_lakesamm + view_otherwater)
# remove non-predictive or unique identifiers
sales_pred <- sales_pred %>%
select(-sale_id, -pinx, -sale_date, -sale_warning, -zoning, -join_status, -join_year, -longitude, -latitude, -city, -area, -land_val, -imp_val, -sqft_1, -sqft_fbsmt, -submarket, -view_sound, -view_lakewash, -view_lakesamm, -view_otherwater, -view_rainier, -view_olympics, -view_cascades, -view_territorial, -view_other, -subdivision, -predicted, -predicted_exp)
# remove non-predictive or unique identifiers
sales_pred <- sales_pred %>%
dplyr::select(-sale_id, -pinx, -sale_date, -sale_warning, -zoning, -join_status, -join_year, -longitude, -latitude, -city, -area, -land_val, -imp_val, -sqft_1, -sqft_fbsmt, -submarket, -view_sound, -view_lakewash, -view_lakesamm, -view_otherwater, -view_rainier, -view_olympics, -view_cascades, -view_territorial, -view_other, -subdivision, -predicted, -predicted_exp)
# change NA to 0 for column sale_nbr
sales_pred <- sales_pred %>%
mutate(sale_nbr = replace_na(sale_nbr, 0))
set.seed(123)
split <- initial_split(sales_pred, prop = 0.8) # 80% training data
train_data <- training(split)
test_data  <- testing(split)
# recipe represent the steps to do data preparation under tidymodels framework
rec <- recipe(sale_price ~ ., data = train_data) %>%
step_normalize(all_numeric_predictors())
lm_spec <- linear_reg() %>%
set_engine("lm")
ridge_spec <- linear_reg(
penalty = tune(),   # 惩罚系数 λ 待调参
mixture = 0         # mixture = 0 表示 ridge
) %>%
set_engine("glmnet")
lasso_spec <- linear_reg(
penalty = tune(),
mixture = 1         # mixture = 1 表示 lasso
) %>%
set_engine("glmnet")
# 4. workflow 集合 ----------------------------------------------------
models <- workflow_set(
preproc = list(base = rec),
models = list(
linear = lm_spec,
ridge  = ridge_spec,
lasso  = lasso_spec
)
)
# 5. 交叉验证设置 -----------------------------------------------------
folds <- vfold_cv(train_data, v = 5)
# 6. 模型拟合 & 调参 --------------------------------------------------
results <- workflow_map(
models,
resamples = folds,
grid = 20,            # 在 20 个 λ 值上交叉验证
metrics = metric_set(rmse, rsq)
)
results
# 7. 查看结果 ---------------------------------------------------------
autoplot(
results,
rank_metric = "rmse",
metric = "rmse"
) +
ggtitle("Comparison of Linear, Ridge, and Lasso Regression") +
facet_wrap(~ wflow_id)  # 🔥 分面显示每种模型
# 8. 汇总表格（平均性能） -------------------------------------------
collect_metrics(results) %>%
select(wflow_id, mean, std_err, .metric) %>%
pivot_wider(names_from = .metric, values_from = mean)
# 8. 汇总表格（平均性能） -------------------------------------------
collect_metrics(results) %>%
select(wflow_id, mean, std_err, .metric) %>%
pivot_wider(names_from = .metric, values_from = mean)
# 8. 汇总表格（平均性能） -------------------------------------------
collect_metrics(results) %>%
dplyr::select(wflow_id, mean, std_err, .metric) %>%
pivot_wider(names_from = .metric, values_from = mean)
# 9. 在测试集上最终评估 ----------------------------------------------
best_lasso <- extract_workflow_set_result(results, "lasso") %>%
select_best(metric = "rmse")
# 8. 汇总表格（平均性能） -------------------------------------------
collect_metrics(results) %>%
dplyr::select(wflow_id, mean, std_err, .metric) %>%
pivot_wider(names_from = .metric, values_from = mean)
# 9. 在测试集上最终评估 ----------------------------------------------
best_lasso <- extract_workflow_set_result(results, "lasso") %>%
dplyr::select_best(metric = "rmse")
# 8. 汇总表格（平均性能） -------------------------------------------
collect_metrics(results) %>%
dplyr::select(wflow_id, mean, std_err, .metric) %>%
pivot_wider(names_from = .metric, values_from = mean)
# 9. 在测试集上最终评估 ----------------------------------------------
best_lasso <- extract_workflow_set_result(results, "lasso") %>%
select_best(metric = "rmse")
sales_pred <- sales %>%
mutate(view_mountains = view_rainier + view_olympics + view_cascades + view_territorial,
view_water = view_sound + view_lakewash + view_lakesamm + view_otherwater)
sales_pred <- sales %>%
mutate(view_mountains = view_rainier + view_olympics + view_cascades + view_territorial,
view_water = view_sound + view_lakewash + view_lakesamm + view_otherwater)
# remove non-predictive features (subjective here)
sales_pred <- sales_pred %>%
dplyr::select(-sale_id, -pinx, -sale_date, -sale_warning, -zoning, -join_status, -join_year, -longitude, -latitude, -city, -area, -land_val, -imp_val, -sqft_1, -sqft_fbsmt, -submarket, -view_sound, -view_lakewash, -view_lakesamm, -view_otherwater, -view_rainier, -view_olympics, -view_cascades, -view_territorial, -view_other, -subdivision, -predicted, -predicted_exp)
# change NA to 0 for column sale_nbr
sales_pred <- sales_pred %>%
mutate(sale_nbr = replace_na(sale_nbr, 0))
# recipe represent the steps to do data preparation under tidymodels framework
rec <- recipe(sale_price ~ ., data = train_data) %>%
update_role(episode_name, new_role = "sale_id") %>%
step_normalize(all_numeric_predictors())
set.seed(123)
split <- initial_split(sales_pred, prop = 0.8) # 80% training data
train_data <- training(split)
test_data  <- testing(split)
# recipe represent the steps to do data preparation under tidymodels framework
rec <- recipe(sale_price ~ ., data = train_data) %>%
update_role(episode_name, new_role = "sale_id") %>%
step_normalize(all_numeric_predictors())
View(train_data)
sales_pred <- sales %>%
mutate(view_mountains = view_rainier + view_olympics + view_cascades + view_territorial,
view_water = view_sound + view_lakewash + view_lakesamm + view_otherwater)
# remove non-predictive features (subjective here)
sales_pred <- sales_pred %>%
dplyr::select(-sale_id, -pinx, -sale_date, -sale_warning, -zoning, -join_status, -join_year, -longitude, -latitude, -city, -area, -land_val, -imp_val, -sqft_1, -sqft_fbsmt, -submarket, -view_sound, -view_lakewash, -view_lakesamm, -view_otherwater, -view_rainier, -view_olympics, -view_cascades, -view_territorial, -view_other, -subdivision, -predicted, -predicted_exp)
# change NA to 0 for column sale_nbr
sales_pred <- sales_pred %>%
mutate(sale_nbr = replace_na(sale_nbr, 0))
set.seed(123)
split <- initial_split(sales_pred, prop = 0.8) # 80% training data
train_data <- training(split)
test_data  <- testing(split)
# recipe represent the steps to do data preparation under tidymodels framework
rec <- recipe(sale_price ~ ., data = train_data) %>%
update_role(episode_name, new_role = "sale_id") %>%
step_normalize(all_numeric_predictors())
View(sales)
sales_pred <- sales %>%
mutate(view_mountains = view_rainier + view_olympics + view_cascades + view_territorial,
view_water = view_sound + view_lakewash + view_lakesamm + view_otherwater)
# remove non-predictive features (subjective here)
sales_pred <- sales_pred %>%
dplyr::select(-pinx, -sale_date, -sale_warning, -zoning, -join_status, -join_year, -longitude, -latitude, -city, -area, -land_val, -imp_val, -sqft_1, -sqft_fbsmt, -submarket, -view_sound, -view_lakewash, -view_lakesamm, -view_otherwater, -view_rainier, -view_olympics, -view_cascades, -view_territorial, -view_other, -subdivision, -predicted, -predicted_exp)
# change NA to 0 for column sale_nbr
sales_pred <- sales_pred %>%
mutate(sale_nbr = replace_na(sale_nbr, 0))
set.seed(123)
split <- initial_split(sales_pred, prop = 0.8) # 80% training data
train_data <- training(split)
test_data  <- testing(split)
# recipe represent the steps to do data preparation under tidymodels framework
rec <- recipe(sale_price ~ ., data = train_data) %>%
update_role(episode_name, new_role = "sale_id") %>%
step_normalize(all_numeric_predictors())
View(train_data)
?step_zv
# recipe represent the steps to do data preparation under tidymodels framework
rec <- recipe(sale_price ~ ., data = train_data) %>%
step_zv(all_numeric(), -all_outcomes()) %>% # remove variables that contain only a single value
step_normalize(all_numeric(), -all_outcomes())
?recipe
rec <- recipe(sale_price ~ ., data = train_data) %>% # model and dataset
step_zv(all_numeric(), -all_outcomes()) %>% # remove variables that contain only a single value
step_normalize(all_numeric(), -all_outcomes()) # normalize (center and scale) the numeric variables
View(rec)
folds <- vfold_cv(train_data, v = 5)
results <- workflow_map(
models,
resamples = folds,
grid = 20,            # 在 20 个 λ 值上交叉验证
metrics = metric_set(rmse, rsq)
)
results
lm_spec <- linear_reg() %>%
set_engine("lm")
lasso_spec <- linear_reg(
penalty = tune(),   # penalty hyperparameter: to be determined later
mixture = 1         # mixture = 1 is lasso and 0 is ridge
) %>%
set_engine("glmnet")
# use `workflow_set()` to put recipt and models together. We have not train yet!
models <- workflow_set(
preproc = list(base = rec),
models = list(
linear = lm_spec,
lasso  = lasso_spec ))
folds <- vfold_cv(train_data, v = 5)
results <- workflow_map(
models,
resamples = folds,
grid = 20,            # 在 20 个 λ 值上交叉验证
metrics = metric_set(rmse, rsq)
)
