Przygotowanie

Załadowanie potrzebnych bibliotek.

library(kableExtra) 
library(reshape2)
library(GGally)
library(ggpmisc)
library(tidyverse)
library(ggpubr)

kable1 = function(data, digits=3) data %>% kable(digits = digits) %>% 
  kable_styling(bootstrap_options = 
                  c("striped", "hover", "condensed", "responsive"))

kable2 = function(data, digits=3, height="320px") data %>% 
  kable1(digits = digits) %>% 
  scroll_box(width = "100%", height = height)

Zadanie 1

Generowanie danych

Dane wygenerowałem do jednej tabeli typu tibble. Pewna wątpliwość pojawiła się jedynie w przypadku zmiennej w7. Zgodnie z polecenie zmienna ta powinna zawierać wartości jak w2, tylko order(w2). Zrozumiałem to w taki sposób, że powinien to być nowy wektor z rozkładu normalnego uporządkowany według porządku z wektora w2. Gdyby bowiem założyć, że chodziło o w7=w2[order(w2)] to tak otrzymany wektor nie różnił by się niczym od sort(w2).

set.seed(1142308)
n=30

df = tibble(
  w1 = runif(n, 1, 20) %>% round(1),
  w2 = rnorm(n, 10, 5) %>% round(2),
  w3 = w1 %>% sort(),
  w4 = w1 %>% sort(decreasing = T),
  w5 = w1 %>% rev(),
  w6 = sample(5:15, n, replace = T),
  w7 = rnorm(n, 10, 5) %>% round(2) %>% .[order(w2)]
)

df %>% kable2()

w1	w2	w3	w4	w5	w6	w7
14.1	12.85	1.4	19.9	7.6	9	15.12
1.5	9.71	1.5	19.4	2.7	8	12.07
16.2	16.15	1.8	18.2	18.2	11	15.15
4.1	10.09	2.7	18.0	3.9	5	6.97
3.8	10.20	3.4	16.2	1.8	14	6.53
10.8	6.41	3.8	15.7	11.2	13	4.88
14.3	7.55	3.9	14.3	13.0	10	4.99
9.7	14.26	4.1	14.1	1.4	5	6.20
6.8	27.61	5.6	13.6	5.6	9	11.79
10.2	11.90	6.2	13.5	13.6	6	6.48
19.9	7.26	6.3	13.0	15.7	14	12.78
19.4	6.92	6.8	11.2	9.5	14	16.01
10.3	5.37	7.6	10.8	3.4	5	16.02
13.5	11.22	9.5	10.3	6.3	13	6.24
18.0	7.24	9.7	10.2	6.2	9	9.99
6.2	2.50	10.2	9.7	18.0	15	12.09
6.3	9.30	10.3	9.5	13.5	6	20.47
3.4	7.38	10.8	7.6	10.3	8	8.28
9.5	14.86	11.2	6.8	19.4	13	8.40
15.7	9.28	13.0	6.3	19.9	5	13.69
13.6	15.34	13.5	6.2	10.2	14	8.02
5.6	13.37	13.6	5.6	6.8	5	14.17
1.4	5.48	14.1	4.1	9.7	14	6.26
13.0	5.71	14.3	3.9	14.3	11	12.91
11.2	9.39	15.7	3.8	10.8	13	3.50
1.8	-1.31	16.2	3.4	3.8	13	11.56
3.9	10.86	18.0	2.7	4.1	10	20.24
18.2	15.79	18.2	1.8	16.2	9	4.60
2.7	14.73	19.4	1.5	1.5	7	12.17
7.6	14.90	19.9	1.4	14.1	13	12.50

Zbadanie korelacji

Korelację przedstawię najpierw w formie przejrzystej mapy cieplnej uzyskanej za pomocą funkcji ggcorr z pakietu GGally. Jako metodę wyznaczenia współczynników korelacji wybrałem metodę rang Spearmana. Wyboru tego dokonałem świadomie mając na uwadze arbitralny brak zgodności z rozkładem normalnym wszystkich zmiennych poza zmienną w2 oraz w7.

df %>% ggcorr(method = c("pairwise", "spearman"), nbreaks = 10, label_round = 2,
              hjust = 0.8, label = TRUE, label_size = 3, color = "grey50")

Największa korekcja jest pomiędzy parami w3-w4, co biorąc pod uwagę metodę generowania tych zmiennych nie jest żadnym zaskoczeniem, oraz w parze w1-w5 co również nie budzi zdziwienia. Jednak te pary nie będą, zgodnie z kolejnymi poleceniami zadania, brane pod uwagę.

Dla wyznaczenia istotności korelacji przygotowałem własna funkcję fcor zwracającą wyznaczone wartości w wygodnej formie tabeli tibble. Niestety w tym przypadku musiałem zrezygnować z metody Spearmana i powrócić do metody Pearsona. W przypadku bowiem metody Spearmana funkcja cor.test nie potrafiła wyznaczyć dokładnej wartości p-value.

fcor = function(df, vars=c(1,2), method="pearson", conf.level = 0.95) {
  if(vars[1]<1 | vars[1]>length(df)) stop("Wybierz własciwy numer kolumny 1")
  if(vars[2]<1 | vars[2]>length(df)) stop("Wybierz własciwy numer kolumny 2")
  
  cortest = cor.test(df[,vars[1]] %>% pull(1), 
                     df[,vars[2]] %>% pull(1), 
                     method=method, conf.level=conf.level)
  tibble(
    pair = paste0(names(df)[vars[1]], "-", names(df)[vars[2]]),
    cor = cortest$estimate,
    t = cortest$statistic,
    df = cortest$parameter,
    p.value = cortest$p.value,
    low = cortest$conf.int[1],
    up = cortest$conf.int[2]
  )
}

Użycie fcor jest bardzo proste.

dfcor = fcor(df, c(1, 2)) %>% 
  bind_rows(fcor(df, c(1, 4))) %>% 
  bind_rows(fcor(df, c(1, 7))) %>% 
  bind_rows(fcor(df, c(3, 6))) %>% 
  bind_rows(fcor(df, c(3, 7))) %>% 
  bind_rows(fcor(df, c(5, 6))) %>% 
  bind_rows(fcor(df, c(2, 7))) 

dfcor %>% kable1()

pair	cor	t	df	p.value	low	up
w1-w2	0.077	0.410	28	0.685	-0.291	0.426
w1-w4	0.131	0.697	28	0.492	-0.241	0.469
w1-w7	-0.034	-0.182	28	0.857	-0.390	0.330
w3-w6	0.131	0.698	28	0.491	-0.241	0.469
w3-w7	0.087	0.465	28	0.646	-0.282	0.434
w5-w6	0.258	1.411	28	0.169	-0.113	0.565
w2-w7	-0.025	-0.131	28	0.896	-0.382	0.338

Spośród wskazanych siedmiu par najwyższy współczynnik korelacji występuje w parze w5-w6 zaś najniższy w przypadku pary w2-w7.

Dodatkowo można zauważyć, iż na poziomie istotności \(\alpha=0.05\) dla każdej z badanych par brak jest podstaw aby odrzucić hipotezę zerowa mówiącą o niezależności zmiennych. Należy więc przyjąć, że w żadnym z badanych przypadków nie występuje istotna korelacja.

Prosta regresji

Zgodnie z poleceniem zadania, wyznaczę prostą regresji dla pary o największym współczynniku korelacji czyli w5-w6 oraz dla pary o najniższym współczynniku korelacji czyli dla pary w2-w7.

Para w5-w6

model = lm(w6~w5, df)
summary(model)

## 
## Call:
## lm(formula = w6 ~ w5, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.5998 -2.0950 -0.3643  2.9880  5.1955 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.5266     1.2319   6.921 1.59e-07 ***
## w5            0.1544     0.1095   1.411    0.169    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.362 on 28 degrees of freedom
## Multiple R-squared:  0.06636,    Adjusted R-squared:  0.03301 
## F-statistic:  1.99 on 1 and 28 DF,  p-value: 0.1694

plot = df %>% ggplot(aes(w5, w6))+
  geom_point(size=3, shape=21, fill="blue")+
  geom_smooth(formula = y~x, method = "lm")+
  stat_correlation(aes(
    label = paste(
      after_stat(cor.label),
      after_stat(t.value.label),
      after_stat(p.value.label),
      after_stat(n.label),
      sep = '*"; "*')), label.y = 0.85, geom = "label_npc")+
  stat_poly_eq(aes(
    label =  paste(
      after_stat(eq.label), "*\" z \"*", 
      after_stat(rr.label), "*\", \"*", 
      after_stat(f.value.label), "*\", oraz \"*",
      after_stat(p.value.label), "*\".\"",
      sep = "")),
    formula = y~x, geom = "label_npc")+
  ylim(c(NA,18))+
  labs(title = "Prosta regresji dla pary zmiennych w5-w6",
       caption = "M. Fiołka")
plot

Dopasowanie modelu jest bardzo słabe. Model wyjaśnia zmienność zmiennej objaśnianej w zaledwie 7 procentach o czym mówi nam współczynnik determinacji \(R^2\). Takie same wnioski należy wyciągnąć z wartości statystyki \(F\) oraz jej p-value która nie pozwala mi na odrzucenie hipotezy zerowej o nieistotności współczynników regresji. Rzeczywiście jak można zauważyć także w przypadku testu \(t\) dla zmiennej w5 również nie ma podstaw aby odrzucić hipotezę zerową mówiącą o zerowej wartości tego współczynnika. Oczywiście, należało się w pełni spodziewać takich wyników już po samej analizie wartości współczynnika korelacji oraz wyniku testu jego istotności.

Para w2-w7

model2 = lm(w7~w2, df)
summary(model2)

## 
## Call:
## lm(formula = w7 ~ w2, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.1912 -4.1540  0.9453  2.7823  9.7769 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.89247    1.89840   5.738 3.72e-06 ***
## w2          -0.02143    0.16306  -0.131    0.896    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.654 on 28 degrees of freedom
## Multiple R-squared:  0.0006166,  Adjusted R-squared:  -0.03508 
## F-statistic: 0.01728 on 1 and 28 DF,  p-value: 0.8964

plot2 = df %>% ggplot(aes(w2, w7))+
  geom_point(size=3, shape=21, fill="blue")+
  geom_smooth(formula = y~x, method = "lm")+
  stat_correlation(aes(
    label = paste(
      after_stat(cor.label),
      after_stat(t.value.label),
      after_stat(p.value.label),
      after_stat(n.label),
      sep = '*"; "*')), label.y = 0.85, geom = "label_npc")+
  stat_poly_eq(aes(
    label =  paste(
      after_stat(eq.label), "*\" z \"*", 
      after_stat(rr.label), "*\", \"*", 
      after_stat(f.value.label), "*\", oraz \"*",
      after_stat(p.value.label), "*\".\"",
      sep = "")),
    formula = y~x, geom = "label_npc")+
  ylim(c(NA,20))+
  labs(title = "Prosta regresji dla pary zmiennych w2-w7",
       caption = "M. Fiołka")
plot2

Próbując wyznaczyć prostą regresji dla pary zmiennych dla których współczynnik korelacji wynosi tylko -0.025 jesteśmy z oczywistych względów skazani na porażkę. Nie dziwi więc ani dramatycznie niska wartość współczynnika determinacji \(R^2\) ani taż bardzo wysokie wartości p-value zarówno w przypadku testu \(F\) Snedecora jak i w przypadku testu \(t\) dla zmiennej objaśniającej, które w tym przypadku miają prawie dokładnie taką samą wartość równą 0.896.

Predykcja

Zbiór predykcyjny (dla pierwszego modelu) będzie przygotowany w formie tabeli danych od razu z standardowym błędem predykcji średniej se.fit, skalkulowanym przedziałem ufności CI (lwr.conf, upr.conf) oraz przedziałem predykcji PI (lwr.pred, upr.pred). Zmienna w6 jest wartością przewidywaną, czyli tym samym co fit i umieszczona została w tej tablicy jedynie w celu późniejszego łatwego dodania przewidywanych wartości na uprzednio stworzonym wykresie.

dfpred = tibble(w5 = sample(1:20, 5))
conf = model %>% predict(dfpred, se.fit=T, interval="confidence", level = 0.95) 
pred = model %>% predict(dfpred, se.fit=T, interval="prediction", level = 0.95)
  
dfpred = dfpred %>% mutate(
  w6 = conf$fit[,1],
  fit = conf$fit[,1],
  se.fit = conf$se.fit,
  lwr.conf = conf$fit[,2],
  upr.conf = conf$fit[,3],
  lwr.pred = pred$fit[,2],
  upr.pred = pred$fit[,3]
) 

dfpred %>% round(1) %>% kable1()

w5	w6	fit	se.fit	lwr.conf	upr.conf	lwr.pred	upr.pred
3	9.0	9.0	1.0	7.0	11.0	1.8	16.2
5	9.3	9.3	0.8	7.6	10.9	2.2	16.4
6	9.5	9.5	0.7	7.9	11.0	2.4	16.5
1	8.7	8.7	1.1	6.3	11.0	1.4	16.0
4	9.1	9.1	0.9	7.3	10.9	2.0	16.3

Teraz łatwo jest nanieść te dane na wcześniej przygotowany wykres.

plot +
   geom_ribbon(aes(ymin = lwr.conf, ymax = upr.conf), 
               dfpred, fill="red", alpha = .2)+ 
   geom_ribbon(aes(ymin = lwr.pred, ymax = upr.pred), 
               dfpred, fill="red", alpha = .1)+ 
  geom_point(data=dfpred, size=3, shape=21, fill="red")+
   ylim(c(NA,21))

Warto zwrócić tutaj uwagę, że przedziały ufności CI (ciemniejsze czerwone wypełnienie) dla tak prognozowanych wartości pokrywają się dokładnie z przedziałem ufności prostej regresji. Natomiast przedziały predykcji PI (czerwone jaśniejsze wypełnienie), uwzględniające wariancję resztową, są już znacznie szersze. W naszym przypadku przedziały te są nawet bardzo szerokie. Jest to niewątpliwie spowodowane bardzo słabym dopasowaniem modelu a co za tym idzie znaczną wariancją rezyduów które zawierają aż 93% zmienności całego modelu. Wszystkie prognozowane punkty leżą oczywiście na naszej prostej regresji, choć w przypadku pierwszego punktu wykonana została nieznaczna i niezamierzona ekstrapolacja.

Zadanie 2

Do przeprowadzenia drugiej części zadania posłużę się danymi o ilości gospodarstw wiejskich oraz powierzchni użytków rolnych w latach 2006-2018 w Polsce. Dane pochodzą z GUS.

Załadowanie danych

library(readxl)

lata = 2006:2018
ROLN_3179 = read_delim(
    "ROLN_3179.csv",
    delim = ";", 
    escape_double = FALSE, 
    trim_ws = TRUE,
    skip = 1,
    col_select = -col29,
    col_names = c(
      "Kod", 
      "Województwo", 
      paste0("gosp_", lata),
      paste0("pow_", lata), 
      "col29"),
    show_col_types = FALSE
    )

ROLN_3179 %>% kable2()

Kod	Województwo	gosp_2006	gosp_2007	gosp_2008	gosp_2009	gosp_2010	gosp_2011	gosp_2012	gosp_2013	gosp_2014	gosp_2015	gosp_2016	gosp_2017	gosp_2018	pow_2006	pow_2007	pow_2008	pow_2009	pow_2010	pow_2011	pow_2012	pow_2013	pow_2014	pow_2015	pow_2016	pow_2017	pow_2018
0000000	POLSKA	9187	11870	14896	17091	20582	23449	25944	26598	24829	22277	22435	20257	19207	228038	287529	314848	367062	519069	605520	661688	669970	657902	580730	536579	494978	484677
0200000	DOLNOŚLĄSKIE	481	652	879	1021	1227	1322	1312	1189	1046	849	813	741	713	19332	21988	28467	26427	39703	45547	44304	37455	37005	31261	29200	27542	27357
0400000	KUJAWSKO-POMORSKIE	173	217	258	279	327	371	390	415	401	363	470	419	395	4846	5884	5943	6826	7688	8376	8813	11152	11573	10645	9263	8331	7655
0600000	LUBELSKIE	1072	1402	1566	1710	1962	2065	2174	2129	1975	1825	1980	1904	1948	19957	23934	26892	30113	34855	34837	37466	40819	38467	34052	31343	29001	28428
0800000	LUBUSKIE	256	361	480	579	833	1081	1356	1422	1370	1202	1148	948	877	12094	18201	18207	22929	35797	44259	52581	54692	53300	46343	43235	37923	37174
1000000	ŁÓDZKIE	218	261	314	366	420	478	518	528	508	478	497	477	491	3378	3561	4829	5270	7671	8746	9908	10342	11229	10158	9986	9260	8905
1200000	MAŁOPOLSKIE	1363	1627	2100	2197	2156	2138	2103	1838	1378	1128	1093	934	770	13827	14481	22619	16488	21968	21396	21050	17005	15529	12976	12364	10691	8844
1400000	MAZOWIECKIE	1028	1215	1481	1673	1935	2140	2373	2609	2374	2147	2426	2215	2284	20878	23218	27742	34538	46229	50100	55804	63445	60354	53790	49517	44348	42049
1600000	OPOLSKIE	46	53	62	63	79	86	90	88	75	67	68	57	61	1196	934	1571	1879	3180	2703	2930	3543	3306	3042	3216	2790	3554
1800000	PODKARPACKIE	1164	1577	1892	2014	2091	2045	1940	1750	1475	1261	1252	1194	1131	20601	27047	28671	22593	31868	32359	30381	29506	23510	16656	15486	15349	13630
2000000	PODLASKIE	628	847	1160	1528	2033	2440	2924	3407	3432	3273	3437	3211	2989	11657	15391	20410	28764	42917	52066	56367	63548	64897	56528	55168	53551	51608
2200000	POMORSKIE	222	273	392	494	648	763	894	893	847	737	679	609	540	8037	10968	11366	14591	22554	27357	30615	28721	29282	24866	23328	22419	19974
2400000	ŚLĄSKIE	116	143	176	199	228	238	236	242	230	201	180	162	148	2340	3079	3935	3632	5739	6787	7125	7220	7788	6638	5324	3726	2951
2600000	ŚWIĘTOKRZYSKIE	892	995	1165	1170	1243	1296	1288	1207	992	853	834	740	680	8963	9824	10841	10647	13123	14301	14551	15122	13038	11598	10739	9970	9087
2800000	WARMIŃSKO-MAZURSKIE	586	773	1059	1514	2279	3033	3793	4235	4234	4041	4142	3745	3393	23991	28810	28791	49617	75242	98473	112945	116199	117097	112768	108667	107067	104574
3000000	WIELKOPOLSKIE	264	415	516	588	748	888	974	1006	966	809	843	736	727	14511	21096	20417	23770	32513	38434	41479	41617	42071	34523	29171	25389	25994
3200000	ZACHODNIOPOMORSKIE	678	1059	1396	1696	2373	3065	3579	3640	3526	3043	2573	2165	2060	42431	59114	54151	68977	98023	119780	135367	129585	129456	114887	100570	87620	92892

Dane pobrane z GUS są w formie zagregowanej i zawierają informację dotyczącą obszaru całego kraju jak i poszczególnych województw. Do dalszej analizy skupię się jedynie na danych dla całego kraju.

dfRoln = ROLN_3179 %>% 
  filter(Województwo == 'POLSKA') %>% 
  select(-Kod, -Województwo) %>% 
  pivot_longer(gosp_2006:pow_2018, names_to = c("name", "rok"),
               names_pattern = "(^[a-z]*)_(\\d*)") %>% 
  pivot_wider() %>% 
  rename(gospodarstwa = gosp,
         powierzchnia = pow) %>% 
  mutate(rok = rok %>% as.integer)

dfRoln %>% kable2()

rok	gospodarstwa	powierzchnia
2006	9187	228038
2007	11870	287529
2008	14896	314848
2009	17091	367062
2010	20582	519069
2011	23449	605520
2012	25944	661688
2013	26598	669970
2014	24829	657902
2015	22277	580730
2016	22435	536579
2017	20257	494978
2018	19207	484677

W ten sposób uzyskałem 13 rekordów dla których zbuduję dwa modele regresji liniowej.

model.gosp = lm(gospodarstwa~rok, dfRoln)
summary(model.gosp)

## 
## Call:
## lm(formula = gospodarstwa ~ rok, data = dfRoln)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5942.3 -3644.5  -244.7  3183.2  6050.0 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -1742396.9   642237.2  -2.713   0.0202 *
## rok              875.9      319.2   2.744   0.0191 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4306 on 11 degrees of freedom
## Multiple R-squared:  0.4064, Adjusted R-squared:  0.3524 
## F-statistic: 7.529 on 1 and 11 DF,  p-value: 0.01909

model.pow = lm(powierzchnia~rok, dfRoln)
summary(model.pow)

## 
## Call:
## lm(formula = powierzchnia ~ rok, data = dfRoln)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -154900  -83266  -52602  116064  168720 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -48669732   18061703  -2.695   0.0209 *
## rok             24435       8977   2.722   0.0199 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 121100 on 11 degrees of freedom
## Multiple R-squared:  0.4025, Adjusted R-squared:  0.3481 
## F-statistic: 7.409 on 1 and 11 DF,  p-value: 0.01986

plot.gosp = dfRoln %>% ggplot(aes(rok, gospodarstwa))+
  geom_point(size=3, shape=21, fill="blue")+
  geom_smooth(formula = y~x, method = "lm")+
  stat_correlation(aes(
    label = paste(
      after_stat(cor.label),
      after_stat(t.value.label),
      after_stat(p.value.label),
      after_stat(n.label),
      sep = '*"; "*')), label.y = 0.85, geom = "label_npc")+
  stat_poly_eq(aes(
    label =  paste(
      after_stat(eq.label), "*\" z \"*", 
      after_stat(rr.label), "*\", \"*", 
      after_stat(f.value.label), "*\", oraz \"*",
      after_stat(p.value.label), "*\".\"",
      sep = "")),
    formula = y~x, geom = "label_npc")+
  labs(title = "Ilość gospodarstw rolnych w Polsce w latach 2006-2018",
       caption = "M. Fiołka")+
    scale_x_continuous(breaks = seq(2006, 2018, 2))
plot.gosp

plot.pow = dfRoln %>% ggplot(aes(rok, powierzchnia))+
  geom_point(size=3, shape=21, fill="blue")+
  geom_smooth(formula = y~x, method = "lm")+
  stat_correlation(aes(
    label = paste(
      after_stat(cor.label),
      after_stat(t.value.label),
      after_stat(p.value.label),
      after_stat(n.label),
      sep = '*"; "*')), label.y = 0.85, geom = "label_npc")+
  stat_poly_eq(aes(
    label =  paste(
      after_stat(eq.label), "*\" z \"*", 
      after_stat(rr.label), "*\", \"*", 
      after_stat(f.value.label), "*\", oraz \"*",
      after_stat(p.value.label), "*\".\"",
      sep = "")),
    formula = y~x, geom = "label_npc")+
  labs(title = "Średnia wielkość godpodarstw rolnych w Polsce w latach 2006-2018",
       caption = "M. Fiołka")+
    scale_x_continuous(breaks = seq(2006, 2018, 2))
plot.pow

W obu przypadkach mamy do czynienia z bardzo słabą liniowością. Ma to swoje odzwierciedlenie w wartości współczynnika determinacji który zarówno dla ilości gospodarstw jak i dla ich wielkości jest prawie identyczny i wynosi 0,4. Co prawda zarówno statystyka \(F\) oraz jej p-value jak i statystyki \(t\) oraz ich wartości p-value pozwalają mi na poziomie \(\alpha=0.05\) na odrzucenie hipotez zerowych o zerowych współczynnikach regresji, ale już nie mógł bym tego zrobić na poziomie \(\alpha=0.01\). Warto jednak zauważyć, że w obu przypadkach w okresie 2012 - 2013 występuje dość zaskakujące maksimum. To skłoniło mnie do podzielenia danych na dwa okresy. Pierwszy okres do roku 2012 włącznie oraz drugi okres od roku 2012 włącznie i zbudowanie modeli na tak rozdzielonych danych.

rok.przełomowy = 2012
dfRoln1 = dfRoln %>% 
  mutate(okres = ifelse(rok<=rok.przełomowy, "okres 1", "okres 2")) %>% 
  bind_rows(dfRoln %>% filter(rok==rok.przełomowy) %>% mutate(okres = "okres 2"))
dfRoln1 %>% kable2()

rok	gospodarstwa	powierzchnia	okres
2006	9187	228038	okres 1
2007	11870	287529	okres 1
2008	14896	314848	okres 1
2009	17091	367062	okres 1
2010	20582	519069	okres 1
2011	23449	605520	okres 1
2012	25944	661688	okres 1
2013	26598	669970	okres 2
2014	24829	657902	okres 2
2015	22277	580730	okres 2
2016	22435	536579	okres 2
2017	20257	494978	okres 2
2018	19207	484677	okres 2
2012	25944	661688	okres 2

Po odpowiednim przygotowaniu danych mogłem przystąpić do budowania modeli dla poszczególnych okresów.

model.gosp1 = lm(gospodarstwa~rok, dfRoln1 %>% filter(okres=="okres 1"))
summary(model.gosp1)

## 
## Call:
## lm(formula = gospodarstwa ~ rok, data = dfRoln1 %>% filter(okres == 
##     "okres 1"))
## 
## Residuals:
##       1       2       3       4       5       6       7 
##   89.46  -53.07  147.39 -483.14  182.32  223.79 -106.75 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.659e+06  1.020e+05  -55.49 3.59e-08 ***
## rok          2.826e+03  5.076e+01   55.67 3.54e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 268.6 on 5 degrees of freedom
## Multiple R-squared:  0.9984, Adjusted R-squared:  0.9981 
## F-statistic:  3099 on 1 and 5 DF,  p-value: 3.538e-08

model.gosp2 = lm(gospodarstwa~rok, dfRoln1 %>% filter(okres=="okres 2"))
summary(model.gosp2)

## 
## Call:
## lm(formula = gospodarstwa ~ rok, data = dfRoln1 %>% filter(okres == 
##     "okres 2"))
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  999.36  490.61 -801.14  617.11 -300.64  -90.39 -914.89 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2562481.9   304527.5   8.415 0.000389 ***
## rok           -1260.2      151.1  -8.339 0.000406 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 799.7 on 5 degrees of freedom
## Multiple R-squared:  0.9329, Adjusted R-squared:  0.9195 
## F-statistic: 69.54 on 1 and 5 DF,  p-value: 0.0004056

model.pow1 = lm(powierzchnia~rok, dfRoln1 %>% filter(okres=="okres 1"))
summary(model.pow1)

## 
## Call:
## lm(formula = powierzchnia ~ rok, data = dfRoln1 %>% filter(okres == 
##     "okres 1"))
## 
## Residuals:
##      1      2      3      4      5      6      7 
##  31197  14218 -34933 -59189  16349  26330   6028 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -153201477   14099710  -10.87 0.000115 ***
## rok              76470       7018   10.90 0.000113 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 37140 on 5 degrees of freedom
## Multiple R-squared:  0.9596, Adjusted R-squared:  0.9515 
## F-statistic: 118.7 on 1 and 5 DF,  p-value: 0.0001131

model.pow2 = lm(powierzchnia~rok, dfRoln1 %>% filter(okres=="okres 2"))
summary(model.pow2)

## 
## Call:
## lm(formula = powierzchnia ~ rok, data = dfRoln1 %>% filter(okres == 
##     "okres 2"))
## 
## Residuals:
##      1      2      3      4      5      6      7 
##  14585  38315  -3059 -11412 -17215   8281 -29495 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 72716471    9413151   7.725 0.000581 ***
## rok           -35798       4672  -7.663 0.000603 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24720 on 5 degrees of freedom
## Multiple R-squared:  0.9215, Adjusted R-squared:  0.9058 
## F-statistic: 58.72 on 1 and 5 DF,  p-value: 0.0006028

plot.gosp1 = dfRoln1 %>% ggplot(aes(rok, gospodarstwa,fill=okres))+
  geom_vline(xintercept = rok.przełomowy, color = "red", linetype="dotdash")+
  geom_point(size=3, shape=21, color="black")+
  stat_poly_line(aes(colour=okres), formula = y~x)+
  stat_poly_eq(aes(
    label =  paste(
      after_stat(eq.label), "*\" z \"*", 
      after_stat(rr.label), "*\", \"*", 
      after_stat(f.value.label), "*\", oraz \"*",
      after_stat(p.value.label), "*\".\"",
      sep = "")),
    formula = y~x, geom = "label_npc", alpha = 0.2)+
    scale_x_continuous(breaks = seq(2006, 2018, 2))+
  ylim(c(NA, 35000))


plot.gosp1 = ggpar(
  plot.gosp1,
  title = "Ilość gospodarstw rolnych w Polsce w latach 2006-2018",
       caption = "M. Fiołka",
       palette = "jco")
plot.gosp1

plot.pow1 = dfRoln1 %>% ggplot(aes(rok, powierzchnia, fill=okres))+
  geom_vline(xintercept = rok.przełomowy, color = "red", linetype="dotdash")+
  geom_point(size=3, shape=21, color="black")+
  stat_poly_line(aes(colour=okres), formula = y~x)+
  stat_poly_eq(aes(
    label =  paste(
      after_stat(eq.label), "*\" z \"*", 
      after_stat(rr.label), "*\", \"*", 
      after_stat(f.value.label), "*\", oraz \"*",
      after_stat(p.value.label), "*\".\"",
      sep = "")),
    formula = y~x, geom = "label_npc", alpha = 0.2)+
    scale_x_continuous(breaks = seq(2006, 2018, 2))+
  ylim(c(NA, 9e5))

plot.pow1 = ggpar(
  plot.pow1,
  title = "Średnia wielkość godpodarstw rolnych w Polsce w latach 2006-2018",
       caption = "M. Fiołka",
       palette = "jco")
plot.pow1

Po rozdzieleniu danych okazało się, że uzyskałem doskonałe dopasowanie w każdym z modeli przy współczynnikach determinacji \(R^2\) bardzo bliskich jedności oraz tak niskich wartościach p-value dla wszystkich testów, że hipotezy zerowe można by odrzucić już na poziomie istotności \(\alpha=0.001\).

Na koniec zadałem sobie pytanie co takiego wydarzyło się w Polsce, że do roku 2012 można obserwować doskolnale liniowy wzrost zarówno ilości gospodarstw rolnych jak i ich powierzchni, po czym tendencja ta się odwróciła i nastąpił doskonale liniowy spadek w obu przypadkach.

Jestem przekonany, że aby wyjaśnić ten efekt należało by zestawić te dane z okresem intensywnego działania pewnej, co by nie mówić bardzo charyzmatycznej postaci polskiej sceny politycznej, lidera rolniczych związków zawodowych a swego czasu nawet wicepremiera i ministra rolnictwa który to z własnego wybory pożegnał się z życiem pod koniec 2011 roku. Jest to dość przerażające, że tylko jedna osoba może mieć aż tak wielki wpływ na przekonania społeczno-ekonomiczne ogromnej części społeczeństwa. Oczywiście jest to tylko i wyłącznie moja własna hipoteza która z pewnością warta jest dalszych, gruntownych badań społeczno-ekonomicznych.

WPM - Mini Projekt

Marek Fiołka

09 03 2022