Файлы для этого занятия можно скачать здесь.

Линейная регрессия с эффектами взаимодействия (linear regression with interactions)

Oценка модели

Возьмем базу, содержащую данные по ценам на квартиры (обсуждалась ранее на курсе “Анализ политологических данных в пакете STATA”).

library(foreign)

df <- read.dta("C:/Users/AllaT/Desktop/Rclass/data_flats.dta")
attach(df)

View(df)

Допустим, мы хотим посмотреть, каким образом на цену квартиры влияет ее жилая площадь (livesp), расстояние до станции метро (metrdist), шаговая доступность (walk), этаж (floor) и при этом учесть эффект взаимодействия переменных “этаж” и “расстояние до станции метро”.

Комментарий к переменным: walk – в шаговой доступности от метро или нет (бинарная), floor – первый/ последний этаж или нет (бинарная).

m1 <- lm(price ~ livesp + metrdist + walk + floor + floor:metrdist, data = df)

# через : указываются переменные, взаимодействие которых учитываем
# не забываем, что те переменные, которые включены в произведение, должны быть
# включены отдельно

summary(m1)
## 
## Call:
## lm(formula = price ~ livesp + metrdist + walk + floor + floor:metrdist, 
##     data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -107.58  -19.86   -3.10   13.72  339.68 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -92.22620    6.01075 -15.344  < 2e-16 ***
## livesp           4.67436    0.09738  48.002  < 2e-16 ***
## metrdist        -1.42559    0.43871  -3.250  0.00117 ** 
## walk            11.84818    1.64787   7.190 9.06e-13 ***
## floor            8.93148    4.47992   1.994  0.04632 *  
## metrdist:floor  -0.07670    0.49286  -0.156  0.87634    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 34.44 on 2034 degrees of freedom
## Multiple R-squared:  0.5604, Adjusted R-squared:  0.5593 
## F-statistic: 518.6 on 5 and 2034 DF,  p-value: < 2.2e-16

Предельные эффекты

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

\(\hat{price}=-92.2+4.7*livesp-1.4*metrdist+11.8*walk+8.9*floor-0.07*metrdist*floor\)

Оценим предельный эффект расстояния до метро на цену квартиры, то есть посчитаем, насколько изменится цена на квартиру при удалении от метро на чуть-чуть (значения остальных переменных не меняются). Для этого посчитаем частную производную по переменной metrdist:

\(\frac{\partial price}{\partial metrdist}=-1.4-0.07*floor\)

Видно, что эффект расстояния на цену квартиры зависит от переменной floor. Как мы знаем, она бинарная, поэтому:

\(floor=0\); \(\frac{\partial price}{\partial metrdist}=-1.4\)

\(floor=1\); \(\frac{\partial price}{\partial metrdist}=-1.4-0.07=-1.47\)

Эффект расстояния до метро на цену квартиры неодинаковый: если квартира находится на первом этаже (floor=1), то с ростом расстояния до станции метро цена на квартиру падает чуть быстрее, чем в случае с квартирой, расположенной на других этажах.

Единственное, стоит помнить, что в нашем примере коэффициент при \(metrdist*floor\) незначимый. Однако часто бывает наоборот: переменные по отдельности не оказывают значимого влияния на отклик, а их “взаимодействие” оказывает.

Конечно, рассчитывать предельные эффекты можно и автоматически. Главное при этом не забывать, что эффекты считаются в определенной точке – при фиксированных значениях переменных (По предыдущему примеру видно, что эффект расстояния на цену квартиры неодинаков при разных значениях floor). Для расчета предельных эффектов в линейных моделях с эффектами взаимодействия нам понадобится пакет LinRegInteractive.

install.packages('LinRegInteractive')

Допустим, мы опять хотим посмотреть на предельный эффект переменной metrdist (расстояние до метро) на цену квартиры в случае, если квартира расположена на первом этаже.

library(LinRegInteractive)
fxInteractive(m1)

После выполнения команды fxInteractive() R предложит выбрать переменную, эффект для которой нужно рассчитать. После указания номера переменной появятся два окна: окно для выставления параметров (Linear Model) и окно с графиком эффектов. В окне Linear Model необходимо выбрать значения переменных, при которых эффект будет считаться, а также выбрать опцию “marginal effect”, если нужен именно предельный эффект. Как мы уже знаем, предельный эффект “metrdist” зависит только от значения переменной floor, поэтому достаточно выставить floor=1, а остальное оставить, как есть по умолчанию. После выставления параметров следует нажать “Snapshot”.

В окне консоли мы увидим большую выдачу, в конце которой есть таблица “Marginal effects in different groups for selected values of metric covariates”, в которой содержится предельный эффект переменной “metrdist” при выставленных значениях остальных переменных.

Call:
lm(formula = price ~ livesp + metrdist + walk + floor + floor:metrdist, 
    data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-107.58  -19.86   -3.10   13.72  339.68 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    -92.22620    6.01075 -15.344  < 2e-16 ***
livesp           4.67436    0.09738  48.002  < 2e-16 ***
metrdist        -1.42559    0.43871  -3.250  0.00117 ** 
walk            11.84818    1.64787   7.190 9.06e-13 ***
floor            8.93148    4.47992   1.994  0.04632 *  
metrdist:floor  -0.07670    0.49286  -0.156  0.87634    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 34.44 on 2034 degrees of freedom
Multiple R-squared:  0.5604,    Adjusted R-squared:  0.5593 
F-statistic: 518.6 on 5 and 2034 DF,  p-value: < 2.2e-16


Selected values of metric covariates
                livesp metrdist      walk floor
value       46.3372549 8.117157 0.6857843     1
ECDF(value)  0.6598039 0.522549 0.3142157     1

Effects in different groups for selected values of metric covariates
   effect
 129.2331

Marginal effects in different groups for selected values of metric covariates
                    livesp  metrdist       walk    floor
value           46.3372549  8.117157  0.6857843 1.000000
ECDF(value)      0.6598039  0.522549  0.3142157 1.000000
marginal effect  4.6743576 -1.502293 11.8481801 8.308856

Конечно, выставлять значения переменных и параметры можно и без помощи окна Linear Model:

# preselect.var - для чего считаем эффект
# preselect.type - 'effect' или 'marginal'
# initial.values - list() значения остальных переменных

fxInteractive(m1, preselect.var = 'metrdist', preselect.type = "effect", initial.values = list(livesp=mean(livesp), walk = 10, floor = 1))

График, который находится во втором окне, не очень наглядный, поэтому лучше подумать и изобразить что-то более понятное и красивое.

Bизуализация эффектов

Можем построить график, чтобы посмотреть, отличается ли влияние расстояния до станции метро на цену квартиры для квартир, расположенных на первом/последнем этаже и на других этажах:

plot(df$metrdist, df$price, type = 'n') # пустой график, чтобы точки не мешали

# проводим регрессионную прямую для каждого типа квартиры (для каждого значения floor)
# и для этого прогоняем модель на разных подвыборках - subset
# плюс при этом меняем цвет прямой

abline(lm(df$price ~ df$metrdist, subset = df$floor==0), lty = 1, col = "red", lwd = 1)
abline(lm(df$price ~ df$metrdist, subset = df$floor==1), lty = 1, col = "green", lwd = 1)
legend("topright", c("floor=0", "floor=1"), col = c("red", "green"), lty = c(1, 1))

Кроме того, можем построить график непосредственно для предельного эффекта какой-либо переменной. Например, опять для metrdist. Для этого нам потребуется пакет interplot.

install.packages("interplot")
library(interplot)
interplot(m1, var1 = "metrdist", var2 = "floor") # модель, эффект чего смотрим (var1), в зависимости от значений какой переменной (var2)

На горизонтальной оси на этом графике отмечены значения переменной floor, на вертикальной – значения предельного эффекта переменной metrdist. Точками на этом графике обозначены значения предельных эффектов переменной metrdist (при floor=0 и floor=1), линиями – доверительные интервалы для этих эффектов (по умолчанию считаются 95%-ные доверительные интервалы).

График выглядит так потому, что переменная floor бинарная. Рассмотрим теперь на другой случай: рассчитаем и покажем на графике предельный эффект переменной floor на цену квартиры в зависимости от значений переменной metrdist.

interplot(m1, var1 = "floor", var2 = "metrdist")

На горизонтальной оси на этом графике отмечены значения переменной metrdist, на вертикальной – значения предельного эффекта переменной floor. Темно-серым цветом закрашены доверительные интервалы для значения предельного эффекта переменной floor на цену квартиры.

Видим, что предельный эффект расположения квартиры (на первом этаже или нет) в целом сильнее, чем эффект расстояния до метро, который мы видели на предыдущем графике. Кроме того, заметно, что влияние расположения квартиры на ее цену несильно зависит от расстояния до ближайшей станции метро: с увеличением расстояния до метро предельный эффект ее расположения изменяется лишь на чуть-чуть. Если мы вернемся к предыдущему графику, то тоже не обнаружим там существенных изменений: эффект расстояния до метро примерно одинаковый как для квартир, расположенных на первом этаже, так и для квартир на всех остальных этажах. Это не случайно – ведь, если вспомнить выдачу с регрессией, коэффициент при переменной взаимодействия metrdist:floor был незначим!

Напоследок, чтобы не запутаться, что на какой оси отмечено, добавим подписи и заголовок к графику. (Графики с помощью interplot строятся как графики с ggplot2, “по слоям”: сам график плюс дополнительные параметры через +):

interplot(m1, var1 = "floor", var2 = "metrdist") + 
  
  xlab("Distance to the closest metro station") +
  
  ylab("Estimated coefficient for location (floor)") +
  
  ggtitle("Marginal effect of \n location of a flat on its price") # \n - перейти на след.строку

Time series and cross-section data

Чтобы понять, в чем разница между TS (time series), CS (cross-section) и TSCS (time series cross-section) данными, рассмотрим три примера:

  1. Есть данные по экономическим показателям в Москве за 1990-2000 годы – TS

  2. Есть данные по экономическим показателям по регионам России за 2014 год – CS

  3. Есть данные по экономическим показателям по регионам России за 1990-2000 годы – TSCS

Фиксированные и случайные эффекты

У нас есть база, содержащая данные по экономическим показателям штатов США за 17 лет (с 1970 по 1986 год).

load("C:/Users/AllaT/Desktop/Rclass/Produc.Rda")
View(Produc)

Для удобства создадим количественную переменную, отвечающую за штат (номер штата):

Produc$nstate <- as.numeric(Produc$state)
attach(Produc)

Теперь хотим посмотреть, как связаны уровень безработицы (unemp) и валовой продукт по штатам (gsp) и при этом учесть, что характер и сила связи между показателями в разных штатах могут отличаться. Построим график, аналогичный тому, который мы строили ранее для цен на квартиры:

plot(unemp, log(gsp), type = 'n') # пустой график, чтобы точки не мешали

# проводим регрессионную прямую для каждого штата (для каждого значения nstate)
# при этом меняем цвет (цвет тоже зависит от i, col=i)
# lwd = 0.5 - толщина линии

for (i in nstate){
abline(lm(log(gsp) ~ unemp, subset = nstate==i), lty = 1, col = i, lwd = 0.5)
}

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

Если мы хотим учесть различия по группам (в данном случае по штатам) и считаем, что влияние каких-то предикторов на отклик примерно одинаковое в разных группах, то нам пригодится модель с фиксированными эффектами. Если считаем, что влияние предикторов на отклик отличается в разных группах, то пригодится модель со случайными эффектами. Однако возможен еще третий вариант: влияние предикторов в разных группах примерно одинаково/неодинаково, но эти сходства/различия не связаны с делением на группы. Иными словами, никаких особых групповых эффектов нет. В таком случае нужно использовать обычную линейную регрессию (OLS).

Конечно, одних наших представлений, каким может быть влияние факторов по группам, и выводов по графику недостаточно для определения типа модели – необходимы формальные критерии. Но прежде чем что-то тестировать, попробуем построить модели разных типов. Для этого нам потребуется пакет plm.

Комментарий: английские эквиваленты: fixed effects model, random effects model, pooled model (OLS).

install.packages('plm')
library(plm)

Mодель с фиксированными эффектами:

pcap – private capital stock, pc – public capital

# fe - просто название модели
# index - группирующая переменная (переменные)
# если модель с фиксированными эффектами - ничего не прописываем дополнительно,
# они в R по умолчанию

fe <- plm(log(gsp) ~ log(pcap) + log(pc) + unemp, data = Produc, index = c("state")) 
summary(fe)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = log(gsp) ~ log(pcap) + log(pc) + unemp, data = Produc, 
##     index = c("state"))
## 
## Balanced Panel: n=48, T=17, N=816
## 
## Residuals :
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -0.16000 -0.03180 -0.00233  0.03100  0.23500 
## 
## Coefficients :
##             Estimate Std. Error  t-value  Pr(>|t|)    
## log(pcap)  0.2932621  0.0355911   8.2398 7.467e-16 ***
## log(pc)    0.7723226  0.0226390  34.1147 < 2.2e-16 ***
## unemp     -0.0193248  0.0011182 -17.2827 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    18.941
## Residual Sum of Squares: 2.059
## R-Squared:      0.8913
## Adj. R-Squared: 0.83559
## F-statistic: 2090.88 on 3 and 765 DF, p-value: < 2.22e-16

Модель со случайными эффектами

# аналогично, но есть аргумент model="random"

re <- plm(log(gsp) ~ log(pcap) + log(pc) + unemp, data = Produc, index = c("state"),model = 'random')
summary(re)
## Oneway (individual) effect Random Effect Model 
##    (Swamy-Arora's transformation)
## 
## Call:
## plm(formula = log(gsp) ~ log(pcap) + log(pc) + unemp, data = Produc, 
##     model = "random", index = c("state"))
## 
## Balanced Panel: n=48, T=17, N=816
## 
## Effects:
##                    var  std.dev share
## idiosyncratic 0.002691 0.051879 0.105
## individual    0.023041 0.151793 0.895
## theta:  0.9174  
## 
## Residuals :
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -0.19800 -0.03270 -0.00199  0.03250  0.26400 
## 
## Coefficients :
##               Estimate Std. Error  t-value  Pr(>|t|)    
## (Intercept) -0.6358474  0.1471273  -4.3217  1.74e-05 ***
## log(pcap)    0.3741636  0.0290958  12.8597 < 2.2e-16 ***
## log(pc)      0.7247214  0.0212900  34.0404 < 2.2e-16 ***
## unemp       -0.0196263  0.0011013 -17.8213 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    24.612
## Residual Sum of Squares: 2.2808
## R-Squared:      0.90733
## Adj. R-Squared: 0.90288
## F-statistic: 2650.08 on 3 and 812 DF, p-value: < 2.22e-16

Обычная линейная регрессия (OLS)

ols <- lm(log(gsp) ~ log(pcap) + log(pc) + unemp, data = Produc)
summary(ols)
## 
## Call:
## lm(formula = log(gsp) ~ log(pcap) + log(pc) + unemp, data = Produc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.45004 -0.09828  0.00593  0.09437  0.58100 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.306720   0.064943  -4.723 2.74e-06 ***
## log(pcap)    0.756201   0.018214  41.517  < 2e-16 ***
## log(pc)      0.337769   0.018614  18.146  < 2e-16 ***
## unemp       -0.010682   0.002567  -4.162 3.49e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.16 on 812 degrees of freedom
## Multiple R-squared:  0.9755, Adjusted R-squared:  0.9755 
## F-statistic: 1.08e+04 on 3 and 812 DF,  p-value: < 2.2e-16

Теперь надо выяснить, какая из моделей действительно нам подходит.

Tестирование модели

F-test

Сравниваем модель с фиксированными эффектами и обычную линейную модель (OLS). Иными словами, проверяем, есть ли групповые эффекты – значимые различия по группам вообще.

Нулевая гипотеза: обычная линейная модель лучше модели с фиксированными эффектами.

pFtest(fe, ols)
## 
##  F test for individual effects
## 
## data:  log(gsp) ~ log(pcap) + log(pc) + unemp
## F = 148, df1 = 47, df2 = 765, p-value < 2.2e-16
## alternative hypothesis: significant effects

Делаем вывод, что есть значимые групповые эффекты. Проверяем дальше.

Hausman test (тест Хаусмана)

Сравниваем модель с фиксированными эффектами и со случайными эффектами.

Нулевая гипотеза: модель со случайными эффектами лучше модели с фиксированными эффектами.

phtest(fe, re)
## 
##  Hausman Test
## 
## data:  log(gsp) ~ log(pcap) + log(pc) + unemp
## chisq = 89.553, df = 3, p-value < 2.2e-16
## alternative hypothesis: one model is inconsistent

Отвергаем нулевую гипотезу и делаем вывод, что модель с фиксированными эффектами лучше.

Если бы модель была со случайными эффектами, то пригодился бы еще один тест.

Breusch-Pagan Lagrange Multiplier test for random effects (тест Бреуша-Пагана)

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

plmtest(re, type = c("bp"))
## 
##  Lagrange Multiplier Test - (Breusch-Pagan)
## 
## data:  log(gsp) ~ log(pcap) + log(pc) + unemp
## chisq = 4565.7, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects

До этого момента мы учитывали только разбиение по штатам. А как быть с годами? Ведь влияние разных показателей на валовой продукт может отличаться не только по штатам, но и по годам. Например, безработица могла оказывать более сильное отрицательное влияние на валовой продукт штатов в 1970-е годы, а в 1980-х ее влияние стало слабее.

Для наглядности будем работать с первыми 10 штатами.

ds <- subset(Produc, nstate<=10)

С одной стороны, можем построить график, аналогичный тому, который мы строили для разных штатов. С другой стороны, можем построить что-то более интересное.

# install.packages("car")

library(car)

# smooth=TRUE - lowess regression (обсуждали на занятии 6)
# reg.line = FALSE - не добавляем обычную регр. прямую
# legend.coord = "right" - добавляем легенду справа

scatterplot(log(gsp) ~ year|nstate, boxplots = FALSE, smooth = TRUE, reg.line = FALSE, 
            data = ds,legend.coord = "right")

По графику видно, что динамика валового продукта у разных штатов примерно одинаковая.

Теперь построим модель с учетом как временных эффектов, так и эффектов по штатам. В данном случае уже однозначно будет использоваться модель с фиксированными эффектами, так как модели со случайными эффектами не годятся для работы с TSCS - данными.

# в index сначала указывается временной показатель (year)

fe2 <- plm(log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, data = Produc, index = c("year","state"))
summary(fe2)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp, 
##     data = Produc, index = c("year", "state"))
## 
## Balanced Panel: n=17, T=48, N=816
## 
## Residuals :
##     Min.  1st Qu.   Median  3rd Qu.     Max. 
## -0.22400 -0.05830 -0.00135  0.04900  0.35800 
## 
## Coefficients :
##             Estimate Std. Error t-value  Pr(>|t|)    
## log(pcap)  0.1647800  0.0174912  9.4207 < 2.2e-16 ***
## log(pc)    0.3035960  0.0104427 29.0727 < 2.2e-16 ***
## log(emp)   0.5888107  0.0137757 42.7428 < 2.2e-16 ***
## unemp     -0.0060575  0.0017702 -3.4220 0.0006534 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    834.46
## Residual Sum of Squares: 6.0629
## R-Squared:      0.99273
## Adj. R-Squared: 0.96719
## F-statistic: 27156 on 4 and 795 DF, p-value: < 2.22e-16

Также можем проверить наличие временных эффектов (опять нулевая гипотеза гласит, что эффектов нет):

plmtest(fe2, c("time"), type = ("bp"))
## 
##  Lagrange Multiplier Test - time effects (Breusch-Pagan)
## 
## data:  log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp
## chisq = 4135, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects

При анализе пространственно-временных данных часто возникает проблема автокорреляции (связь между остатками, относящихся к наблюдениям в разных временных группах, serial correlation). Кроме того, может возникать другая проблема – связь между остатками, относящихся к наблюдениям в разных “пространственных” группах (cross-sectional dependence). Обычно считается, что эти явления действительно представляют проблему, когда временных групп (лет) больше, чем пространственных (страны, регионы, штаты). В нашем случае есть данные по 48 штатам за 13 лет. Проверим:

pcdtest(fe2, test = c("lm"))
## 
##  Breusch-Pagan LM test for cross-sectional dependence in panels
## 
## data:  formula
## chisq = 4651.3, df = 136, p-value < 2.2e-16
## alternative hypothesis: cross-sectional dependence
pcdtest(fe2, test = c("cd"))
## 
##  Pesaran CD test for cross-sectional dependence in panels
## 
## data:  formula
## z = 67.449, p-value < 2.2e-16
## alternative hypothesis: cross-sectional dependence

Видно, что связь между остатками для наблюдений, относящихся к разным штатам, есть. Это мы тестировали наличие cross-sectional dependence. Теперь посмотрим, связаны ли остатки у наблюдений, относящихся к разным годам (serial correlation):

pbgtest(fe2)
## 
##  Breusch-Godfrey/Wooldridge test for serial correlation in panel
##  models
## 
## data:  log(gsp) ~ log(pcap) + log(pc) + log(emp) + unemp
## chisq = 708.08, df = 48, p-value < 2.2e-16
## alternative hypothesis: serial correlation in idiosyncratic errors

Видим, что автокорреляция присутствует.

Чтобы избежать смещенных оценок коэффициентов, скорректируем оценку модели на наличие гетероскедастичности и автокорреляции – при оценивании модели будем использовать не обычные стандартные ошибки, а устойчивые к гетероскедастичности и автокорреляции (hac standard errors).

Воспользуемся уже известным пакетом lmtest.

# install.packages("lmtest")

library(lmtest)

# "arellano" - учитываем как наличие гетероскедастичности, так и автокорреляции

coeftest(fe2, vcovHC(fe2, method = "arellano", type = "HC3"))
## 
## t test of coefficients:
## 
##             Estimate Std. Error t value Pr(>|t|)    
## log(pcap)  0.1647800  0.0183272   8.991  < 2e-16 ***
## log(pc)    0.3035960  0.0088167  34.434  < 2e-16 ***
## log(emp)   0.5888107  0.0228968  25.716  < 2e-16 ***
## unemp     -0.0060575  0.0029235  -2.072  0.03859 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Получили скорректированные оценки коэффициентов.

Подробнее о работе с пакетом plm для моделей с эффектами можно почитать, например, здесь.