Машинное обучение со смешанными эффектами с GPBoost для сгруппированных и площадных пространственных эконометрических данных

Демонстрация с использованием данных о ВВП Европы

Алгоритм GPBoost расширяет модели линейных смешанных эффектов и гауссовских процессов, заменяя линейную функцию с фиксированными эффектами непараметрической нелинейной функцией, смоделированной с использованием повышения дерева. В данной статье показано, как алгоритм GPBoost, реализованный в библиотекеGPBoost, можно использовать для моделирования данных с пространственной и групповой структурой. Мы демонстрируем функциональность библиотеки GPBoost, используя данные европейского ВВП, которые являются примером площадных пространственных эконометрических данных.

Введение

∘ Введение
· Оглавление
· Основной рабочий процесс GPBoost
· Описание данных
· Загрузка данных и краткая визуализация
∘ Обучение модели GPBoost
∘ Выбор параметров настройки
∘ Интерпретация модели
· Модель с оценкой случайных эффектов
· Пространственная карта эффектов
· Понимание функции фиксированных эффектов
∘ Расширения
· Разделение случайных эффектов на разное время периоды
· Взаимодействие между пространственными и предикторными переменными с фиксированными эффектами
· Большие данные
· Другие модели пространственных случайных эффектов
· (Обобщенные) линейные смешанные эффекты и модели гауссовских процессов
∘Ссылки

Основной рабочий процесс GPBoost

Применение модели GPBoost (= комбинированное повышение дерева и случайные эффекты / модели GP) включает следующие основные этапы:

Определите GPModel, в котором вы указываете следующее:
— модель случайных эффектов (например, пространственные случайные эффекты, сгруппированные случайные эффекты, комбинированные пространственные и сгруппированные и т. д.)
— вероятность (= распределение переменная отклика, зависящая от фиксированных и случайных эффектов)
Создайте gpb.Dataset с переменной ответа и фиксированными переменными-предикторами эффектов
Выберите параметры настройки, например, с помощью функции gpb.grid.search.tune.parameters
Обучите модель функцией gpboost / gpb.train
Интерпретировать обученную модель и/или делать прогнозы

В этой статье мы демонстрируем эти шаги, используя реальный набор данных. Кроме того, мы также покажем, как интерпретировать подобранные модели, и рассмотрим несколько расширений и дополнительных функций библиотеки GPBoost. Все результаты получены с использованием GPBoost версии 1.2.1. В этой демонстрации используется пакет R, но соответствующий пакет Python обеспечивает ту же функциональность.

Описание данных

Данные, используемые в этой демонстрации, представляют собой европейские данные о валовом внутреннем продукте (ВВП). Его можно загрузить с https://raw.githubusercontent.com/fabsig/GPBoost/master/data/gdp_european_regions.csv. Данные были собраны Массимо Джаннини из Римского университета Тор Вергата из Евростата и любезно предоставлены Фабио Зигристу для выступления в Римском университете Тор Вергата 16 июня 2023 года.

Данные были собраны по 242 европейским регионам за два года, 2000 и 2021. То есть общее количество точек данных составляет 484. Переменная ответа — (логарифм) ВВП на душу населения. Есть четыре переменные-предикторы:

L: (log) доля занятости (empl/pop)
K: (log) основной капитал/население
Pop: журнал (население)
Edu: доля высшего образования

Кроме того, имеются пространственные координаты центра тяжести каждого региона (долгота и широта), идентификатор пространственного региона для 242 различных европейских регионов и дополнительный идентификатор пространственного кластера, который идентифицирует кластер, к которому принадлежит регион (существует два кластера).

Загрузка данных и краткая визуализация

Сначала мы загружаем данные и создаем карту, иллюстрирующую (логарифмический) ВВП на душу населения в пространстве. Мы создаем две карты: одну со всеми данными и другую при исключении некоторых удаленных островов. В приведенном ниже коде с комментариями мы также показываем, как создать пространственный график, когда файл формы для пространственных полигонов недоступен.

library(gpboost)
library(ggplot2)
library(gridExtra)
library(viridis)
library(sf)

# Load data
data <- read.csv("https://raw.githubusercontent.com/fabsig/GPBoost/master/data/gdp_european_regions.csv")
FID <- data$FID
data <- as.matrix(data[,names(data)!="FID"]) # convert to matrix since the boosting part currently does not support data.frames
covars <- c("L", "K", "pop", "edu")
# Load shape file for spatial plots
cur_tempfile <- tempfile()
download.file(url = "https://raw.githubusercontent.com/fabsig/GPBoost/master/data/shape_European_regions.zip", destfile = cur_tempfile)
out_directory <- tempfile()
unzip(cur_tempfile, exdir = out_directory)
shape <- st_read(dsn = out_directory)

# Create spatial plot of GDP
data_plot <- merge(shape, data.frame(FID = FID, y = data[,"y"]), by="FID")
p1 <- ggplot(data_plot) + geom_sf(aes(group=FID, fill=y)) +
  scale_fill_viridis(name="GDP / capita (log)", option = "B") + 
   ggtitle("GDP / capita (log)")
# Sample plot excluding islands
p2 <- ggplot(data_plot) + geom_sf(aes(group=FID, fill=y)) +
  scale_fill_viridis(name="GDP / capita (log)", option = "B") +
  xlim(2700000,6400000) + ylim(1500000,5200000) + 
  ggtitle("GDP / capita (log) -- excluding islands")
grid.arrange(p1, p2, ncol=2)

# # Spatial plot without a shape file
# p1 <- ggplot(data = data.frame(Lat=data[,"Lat"], Long=data[,"Long"],
#                                GDP=data[,"y"]), aes(x=Long,y=Lat,color=GDP)) +
#   geom_point(size=2, alpha=0.5) + scale_color_viridis(option = "B") + 
#   ggtitle("GDP / capita (log)")
# p2 <- ggplot(data = data.frame(Lat=data[,"Lat"], Long=data[,"Long"], 
#                                GDP=data[,"y"]), aes(x=Long,y=Lat,color=GDP)) +
#   geom_point(size=3, alpha=0.5) + scale_color_viridis(option = "B") + 
#   ggtitle("GDP / capita (log) -- Europe excluding islands") + xlim(-10,28) + ylim(35,67)
# grid.arrange(p1, p2, ncol=2)

Обучение модели GPBoost

Далее мы используем модель гауссовского процесса с экспоненциальной ковариационной функцией для моделирования пространственных случайных эффектов. Кроме того, мы включаем сгруппированные случайные эффекты для переменной кластера cl. В библиотеке GPBoost случайные эффекты процесса Гаусса определяются аргументом gp_coords, а сгруппированные случайные эффекты — аргументом group_data конструктора GPModel. Вышеупомянутые переменные-предикторы используются в функции ансамбля деревьев с фиксированными эффектами. Мы подгоняем модель, используя функцию gpboost или, что то же самое, функцию gpb.train. Обратите внимание, что мы используем параметры настройки, выбранные ниже с помощью перекрестной проверки.

gp_model <- GPModel(group_data = data[, c("cl")], 
                    gp_coords = data[, c("Long", "Lat")],
                    likelihood = "gaussian", cov_function = "exponential")
params <- list(learning_rate = 0.01, max_depth = 2, num_leaves = 2^10,
               min_data_in_leaf = 10, lambda_l2 = 0)
nrounds <- 37
# gp_model$set_optim_params(params = list(trace=TRUE)) # To monitor hyperparameter estimation
boost_data <- gpb.Dataset(data = data[, covars], label = data[, "y"])
gpboost_model <- gpboost(data = boost_data, gp_model = gp_model, 
                         nrounds = nrounds, params = params, 
                         verbose = 1) # same as gpb.train# same as gpb.train gpb.train

Выбор параметров настройки

Важно, чтобы параметры настройки были правильно выбраны для усиления. Универсальных значений по умолчанию не существует, и для каждого набора данных, вероятно, потребуются разные параметры настройки. Ниже мы покажем, как можно выбрать параметры настройки с помощью функции gpb.grid.search.tune.parameters. Мы используем среднеквадратичную ошибку (mse) в качестве меры точности прогноза для данных проверки. В качестве альтернативы можно также использовать, например, тест отрицательного логарифмического правдоподобия (test_neg_log_likelihood = значение по умолчанию, если ничего не указано), который также учитывает неопределенность прогноза.В зависимости от набора данных и размера сетки это может занять некоторое время. время. Вместо детерминированного поиска в сетке, как показано ниже, можно также выполнить случайный поиск в сетке, чтобы ускорить процесс (см. num_try_random).

gp_model <- GPModel(group_data = data[, "cl"], 
                    gp_coords = data[, c("Long", "Lat")],
                    likelihood = "gaussian", cov_function = "exponential")
boost_data <- gpb.Dataset(data = data[, covars], label = data[, "y"])
param_grid = list("learning_rate" = c(1,0.1,0.01), 
                  "min_data_in_leaf" = c(10,100,1000),
                  "max_depth" = c(1,2,3,5,10),
                  "lambda_l2" = c(0,1,10))
other_params <- list(num_leaves = 2^10)
set.seed(1)
opt_params <- gpb.grid.search.tune.parameters(param_grid = param_grid, params = other_params,
                                              num_try_random = NULL, nfold = 4,
                                              data = boost_data, gp_model = gp_model, 
                                              nrounds = 1000, early_stopping_rounds = 10,
                                              verbose_eval = 1, metric = "mse") # metric = "test_neg_log_likelihood"
opt_params
# ***** New best test score (l2 = 0.0255393919591794) found for the following parameter combination: learning_rate: 0.01, min_data_in_leaf: 10, max_depth: 2, lambda_l2: 0, nrounds: 37

Интерпретация модели

Расчетная модель случайных эффектов

Информацию о расчетной модели случайных эффектов можно получить, вызвав функцию summary для gp_model. Для гауссовских процессов GP_var — это предельная дисперсия, т. е. величина пространственной корреляции или пространственной вариации структуры, а GP_range — это параметр диапазона или масштаба, который измеряет, насколько быстро корреляция затухает в пространстве. Для экспоненциальной ковариационной функции утроенное это значение (здесь примерно 17) представляет собой расстояние, на котором (остаточная) пространственная корреляция практически равна нулю (ниже 5%). Как показывают приведенные ниже результаты, величина пространственной корреляции относительно невелика, поскольку предельная дисперсия 0,0089 мала по сравнению с общей дисперсией переменной отклика, которая составляет прибл. 0,29. Кроме того, член ошибки и сгруппированные случайные эффекты также имеют небольшую дисперсию (0,013 и 0,022). Таким образом, мы заключаем, что большая часть дисперсии переменной отклика объясняется переменными-предикторами с фиксированными эффектами.

summary(gp_model)
## =====================================================
## Covariance parameters (random effects):
##            Param.
## Error_term 0.0131
## Group_1    0.0221
## GP_var     0.0089
## GP_range   5.7502
## =====================================================

Карта пространственного эффекта

Мы можем построить предполагаемые («прогнозированные») пространственные случайные эффекты в местоположениях обучающих данных, вызвав функцию predict для обучающих данных; см. код ниже. Такой график показывает пространственный эффект при факторинге эффекта переменных-предикторов с фиксированными эффектами. Обратите внимание, что эти пространственные эффекты учитывают как пространственный гауссовский процесс, так и случайные эффекты кластера сгруппированных областей. Если кто-то хочет получить только пространственные случайные эффекты от части гауссовского процесса, можно использовать функцию predict_training_data_random_effects (см. код с комментариями ниже). В качестве альтернативы можно также выполнить пространственную экстраполяцию (=«Кригинг»), но это не имеет особого смысла для площадных данных.

pred <- predict(gpboost_model, group_data_pred = data[1:242, c("cl")], 
                gp_coords_pred = data[1:242, c("Long", "Lat")],
                data = data[1:242, covars], predict_var = TRUE, pred_latent = TRUE)
data_plot <- merge(shape, data.frame(FID = FID, y = pred$random_effect_mean), by="FID")
plot_mu <- ggplot(data_plot) + geom_sf(aes(group=FID, fill=y)) +
  scale_fill_viridis(name="Spatial effect", option = "B") +
  ggtitle("Spatial effect (mean)") + xlim(2700000,6400000) + ylim(1500000,5200000)
data_plot <- merge(shape, data.frame(FID = FID, y = pred$random_effect_cov), by="FID")
plot_sd <- ggplot(data_plot) + geom_sf(aes(group=FID, fill=y)) +
  scale_fill_viridis(name="Std. dev.", option = "B") +
  ggtitle("Uncertainty (std. dev.)") + xlim(2700000,6400000) + ylim(1500000,5200000)
grid.arrange(plot_mu, plot_sd, ncol=2)

# Only spatial effecst from the Gaussian process
# rand_effs <- predict_training_data_random_effects(gp_model, predict_var = TRUE)[1:242, c("GP", "GP_var")]
# data_plot <- merge(shape, data.frame(FID = FID, y = rand_effs[,1]), by="FID")
# plot_mu <- ggplot(data_plot) + geom_sf(aes(group=FID, fill=y)) +
#   scale_fill_viridis(name="Spatial effect (mean)", option = "B") +
#   ggtitle("Spatial effect from Gausian process (mean)") + xlim(2700000,6400000) + ylim(1500000,5200000)
# data_plot <- merge(shape, data.frame(FID = FID, y = rand_effs[,2]), by="FID")
# plot_sd <- ggplot(data_plot) + geom_sf(aes(group=FID, fill=y)) +
#   scale_fill_viridis(name="Uncertainty (std. dev.)", option = "B") +
#   ggtitle("Uncertainty (std. dev.)") + xlim(2700000,6400000) + ylim(1500000,5200000)
# grid.arrange(plot_mu, plot_sd, ncol=2)

Понимание функции фиксированных эффектов

Существует несколько инструментов для понимания формы функции фиксированных эффектов. Ниже мы рассмотрим меры переменной важности, меры взаимодействия и графики зависимости. В частности, мы смотрим на

SHAP-значения
Графики зависимости SHAP
Важность переменной на основе разделения
H-статистика Фридмана
Графики частных зависимостей (одномерные и двумерные)
График силы SHAP

Как показывают приведенные ниже результаты, информация, полученная из значений SHAP и графиков зависимости SHAP, такая же, как и при просмотре традиционных показателей важности переменных и графиков частичной зависимости. Наиболее важными переменными являются «K» и «edu». Из графиков зависимостей делаем вывод о наличии нелинейностей. Например, эффект K почти плоский для больших и малых значений K и увеличивается в промежутке между ними. Кроме того, влияние edu становится более резким при малых значениях edu и сглаживается при больших значениях edu. H-статистика Фридмана указывает на наличие некоторых взаимодействий. Для двух переменных с наибольшим количеством взаимодействий, L и pop, мы создаем двумерный график частичной зависимости ниже. Кроме того, график силы SHAP показывает, что влияние переменных-предикторов отличается для 2000 года по сравнению с 2021 годом.

# SHAP values
library(SHAPforxgboost)
shap.plot.summary.wrap1(gpboost_model, X = data[,covars]) + ggtitle("SHAP values")

# SHAP dependence plots
shap_long <- shap.prep(gpboost_model, X_train = data[,covars])
shap.plot.dependence(data_long = shap_long, x = covars[2], 
                           color_feature = covars[4], smooth = FALSE, size = 2) + 
  ggtitle("SHAP dependence plot for K")
shap.plot.dependence(data_long = shap_long, x = covars[4], 
                           color_feature = covars[2], smooth = FALSE, size = 2) + 
  ggtitle("SHAP dependence plot for edu")

# SHAP force plot
shap_contrib <- shap.values(gpboost_model, X_train = data[,covars])
plot_data <- cbind(shap_contrib$shap_score, ID = 1:dim(data)[1])
shap.plot.force_plot(plot_data, zoom_in=FALSE, id = "ID") + 
  scale_fill_discrete(name="Variable") + xlab("Regions") + 
  geom_vline(xintercept=242.5, linewidth=1) + ggtitle("SHAP force plot") + 
  geom_text(x=100,y=-1.2,label="2000") + geom_text(x=400,y=-1.2,label="2021")

# Split-based feature importances
feature_importances <- gpb.importance(gpboost_model, percentage = TRUE)
gpb.plot.importance(feature_importances, top_n = 25, measure = "Gain", 
                    main = "Split-based variable importances")

# H-statistic for interactions
library(flashlight)
fl <- flashlight(model = gpboost_model, data = data.frame(y = data[,"y"], data[,covars]), 
                 y = "y", label = "gpb",
                 predict_fun = function(m, X) predict(m, data.matrix(X[,covars]), 
                                                      gp_coords_pred = matrix(-100, ncol = 2, nrow = dim(X)[1]),
                                                      group_data_pred = matrix(-1, ncol = 1, nrow = dim(X)[1]),
                                                      pred_latent = TRUE)$fixed_effect)
plot(imp <- light_interaction(fl, v = covars, pairwise = TRUE)) + 
  ggtitle("H interaction statistic") # takes a few seconds

# Partial dependence plots
gpb.plot.partial.dependence(gpboost_model, data[,covars], variable = 2, xlab = covars[2], 
                            ylab = "gdp", main = "Partial dependence plot" )
gpb.plot.partial.dependence(gpboost_model, data[,covars], variable = 4, xlab = covars[4], 
                            ylab = "gdp", main = "Partial dependence plot" )

# Two-dimensional partial dependence plot (to visualize interactions)
i = 1; j = 3;# i vs j
gpb.plot.part.dep.interact(gpboost_model, data[,covars], variables = c(i,j), xlab = covars[i], 
                           ylab = covars[j], main = "Pairwise partial dependence plot")

Расширения

Отдельные случайные эффекты для разных периодов времени

В приведенной выше модели мы использовали одни и те же случайные эффекты для 2000 и 2021 годов. В качестве альтернативы можно также использовать независимые пространственные и сгруппированные случайные эффекты для разных периодов времени (независимые при априорных, обусловленных данными, зависимость…). В библиотеке GPBoost это можно сделать с помощью аргумента cluster_ids. cluster_ids должен быть вектором длины, равной размеру выборки, и каждая запись указывает на «кластер», к которому принадлежит наблюдение. Затем разные кластеры имеют независимые пространственные и сгруппированные случайные эффекты, но гиперпараметры (например, предельная дисперсия, компоненты дисперсии и т. д.) и функция фиксированных эффектов одинаковы для всех кластеров.

Ниже мы покажем, как мы можем подогнать такую модель и создать две отдельные пространственные карты. Как показывают приведенные ниже результаты, пространственные эффекты для двух лет, 2000 и 2021 гг., сильно различаются. В частности, для 2021 г. наблюдается меньшая (остаточная) пространственная изменчивость (или неоднородность) по сравнению с 2000 г. Это подтверждается рассмотрением стандартных отклонения предсказанных пространственных случайных эффектов, которые почти вдвое больше для 2000 г. по сравнению с 2021 г. (см. ниже). Еще один вывод состоит в том, что в 2000 г. было больше регионов с «низкими» показателями ВВП (пространственные эффекты ниже 0), а в 2021 г. это уже не так.

gp_model <- GPModel(group_data = data[, c("cl")], gp_coords = data[, c("Long", "Lat")],
                    likelihood = "gaussian", cov_function = "exponential",
                    cluster_ids  = c(rep(1,242), rep(2,242)))
boost_data <- gpb.Dataset(data = data[, covars], label = data[, "y"])
params <- list(learning_rate = 0.01, max_depth = 1, num_leaves = 2^10,
               min_data_in_leaf = 10, lambda_l2 = 1) 
# Note: we use the same tuning parameters as above. Ideally, the would have to be chosen again
gpboost_model <- gpboost(data = boost_data, gp_model = gp_model, nrounds = nrounds,
                         params = params, verbose = 0)
# Separate spatial maps for the years 2000 and 2021
pred <- predict(gpboost_model, group_data_pred = data[, c("cl")], 
                gp_coords_pred = data[, c("Long", "Lat")],
                data = data[, covars], predict_var = TRUE, pred_latent = TRUE,
                cluster_ids_pred = c(rep(1,242), rep(2,242)))
data_plot <- merge(shape, data.frame(FID = FID, y = pred$random_effect_mean[1:242]), by="FID")
plot_mu_2000 <- ggplot(data_plot) + geom_sf(aes(group=FID, fill=y)) +
  scale_fill_viridis(name="Spatial effect", option = "B") +
  ggtitle("Spatial effect for 2000 (mean)") + xlim(2700000,6400000) + ylim(1500000,5200000)
data_plot <- merge(shape, data.frame(FID = FID, y = pred$random_effect_mean[243:484]), by="FID")
plot_mu_2021 <- ggplot(data_plot) + geom_sf(aes(group=FID, fill=y)) +
  scale_fill_viridis(name="Spatial effect", option = "B") +
  ggtitle("Spatial effect for 2021 (mean)") + xlim(2700000,6400000) + ylim(1500000,5200000)
grid.arrange(plot_mu_2000, plot_mu_2021, ncol=2)
sd(pred$random_effect_mean[1:242]) # 0.2321267
sd(pred$random_effect_mean[243:484]) # 0.1286398

Взаимодействие между пространственными и предикторными переменными с фиксированными эффектами

В приведенной выше модели нет взаимодействия между случайными эффектами и предикторными переменными с фиксированными эффектами. То есть между пространственными координатами и предикторными переменными с фиксированными эффектами нет взаимодействия. Такое взаимодействие можно смоделировать, дополнительно включив входные переменные случайных эффектов (= пространственные координаты или категориальную группирующую переменную) в функцию фиксированных эффектов. Код ниже показывает, как это можно сделать. Как показано на приведенном ниже графике важности переменных, координаты не используются в древовидном ансамбле, и поэтому для этого набора данных такое взаимодействие не обнаруживается.

gp_model <- GPModel(group_data = data[, c("cl")], gp_coords = data[, c("Long", "Lat")],
                    likelihood = "gaussian", cov_function = "exponential")
covars_interact <- c(c("Long", "Lat"), covars) ## add spatial coordinates to fixed effects predictor variables
boost_data <- gpb.Dataset(data = data[, covars_interact], label = data[, "y"])
params <- list(learning_rate = 0.01, max_depth = 1, num_leaves = 2^10,
               min_data_in_leaf = 10, lambda_l2 = 1) 
# Note: we use the same tuning parameters as above. Ideally, the would have to be chosen again
gpboost_model <- gpboost(data = boost_data, gp_model = gp_model, nrounds = nrounds,
                         params = params, verbose = 0)
feature_importances <- gpb.importance(gpboost_model, percentage = TRUE)
gpb.plot.importance(feature_importances, top_n = 25, measure = "Gain", 
                    main = "Variable importances when including coordinates in the fixed effects")

Большие данные

Для больших наборов данных вычисления с гауссовскими процессами выполняются медленно, и приходится использовать приближение. Библиотека GPBoost реализует несколько из них. Например, установка gp_approx="vecchia" в конструкторе GPModel будет использовать приближение Веккья. Набор данных, использованный в этой статье, относительно невелик, и мы можем точно выполнить все расчеты.

Другие модели пространственных случайных эффектов

Выше мы использовали гауссовский процесс для моделирования пространственных случайных эффектов. Поскольку данные представляют собой территориальные данные, другим вариантом является использование моделей, основанных на информации о районе, таких как модели CAR и SAR. В настоящее время эти модели еще не реализованы в библиотеке GPBoost (могут быть добавлены в будущем —› свяжитесь с автором).

Другой вариант — использовать модель сгруппированных случайных эффектов, определяемую категориальной переменной идентификатора региона, для моделирования пространственных эффектов. Код ниже показывает, как это можно сделать в GPBoost. Однако эта модель по существу игнорирует пространственную структуру.

gp_model <- GPModel(group_data = data[, c("group", "cl")], likelihood = "gaussian")

(Обобщенные) линейные смешанные эффекты и модели гауссовых процессов

(Обобщенные) линейные смешанные эффекты и модели гауссовых процессов также могут быть запущены в библиотеке GPBoost. В приведенном ниже коде показано, как это можно сделать с помощью функции fitGPModel. Обратите внимание, что нужно вручную добавить столбец из 1, чтобы включить перехват.

X_incl_1 = cbind(Intercept = rep(1,dim(data)[1]), data[, covars])
gp_model <- fitGPModel(group_data = data[, c("cl")], gp_coords = data[, c("Long", "Lat")],
                       likelihood = "gaussian", cov_function = "exponential",
                       y = data[, "y"], X = X_incl_1, params = list(std_dev=TRUE))
summary(gp_model)
## =====================================================
## Model summary:
## Log-lik     AIC     BIC 
##  195.97 -373.94 -336.30 
## Nb. observations: 484
## Nb. groups: 2 (Group_1)
## -----------------------------------------------------
## Covariance parameters (random effects):
##            Param. Std. dev.
## Error_term 0.0215    0.0017
## Group_1    0.0003    0.0008
## GP_var     0.0126    0.0047
## GP_range   7.2823    3.6585
## -----------------------------------------------------
## Linear regression coefficients (fixed effects):
##            Param. Std. dev. z value P(>|z|)
## Intercept 16.2816    0.4559 35.7128  0.0000
## L          0.4243    0.0565  7.5042  0.0000
## K          0.6493    0.0212 30.6755  0.0000
## pop        0.0134    0.0097  1.3812  0.1672
## edu        0.0100    0.0009 10.9645  0.0000
## =====================================================