project1<-read.csv("C:/Users/Allan/OneDrive/Desktop/taluremia clean data analysis 1.csv")
project1<-na.omit(project1)
#selecting data for NE region
NE_region<-project1[project1$county=="NE",]
# modelling for NE region
library(plm)
names(project1)
## [1] "county" "state" "state.1" "year"
## [5] "year.1" "forest" "mixeddev" "CumTemp"
## [9] "CumPrecip" "TempVar" "PrecipVar" "AvgWinterTemp"
## [13] "AvgSpringPrecip" "HotDryDays" "talremia_inc" "tick_lag"
NE_region_model1<-plm(talremia_inc~TempVar+tick_lag,data=NE_region,index=c("county","year.1"),model="within")
## Warning in pdata.frame(data, index = index, ...): duplicate couples (id-time) in resulting pdata.frame
## to find out which, use, e.g., table(index(your_pdataframe), useNA = "ifany")
summary(NE_region_model1)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = talremia_inc ~ TempVar + tick_lag, data = NE_region,
## model = "within", index = c("county", "year.1"))
##
## Balanced Panel: n = 1, T = 10, N = 110
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -2.96859 -1.08948 -0.33076 0.54145 7.68691
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## TempVar -0.0143082 0.0028091 -5.0935 1.519e-06 ***
## tick_lag 0.0289948 0.0178801 1.6216 0.1078
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 419.42
## Residual Sum of Squares: 333.63
## R-Squared: 0.20454
## Adj. R-Squared: 0.18967
## F-statistic: 13.7568 on 2 and 107 DF, p-value: 4.82e-06
NE_region_model2<-plm(talremia_inc ~CumPrecip+forest
+mixeddev
+TempVar
+PrecipVar
+AvgWinterTemp
+AvgSpringPrecip
+HotDryDays
+tick_lag,data=NE_region,index=c("county","year.1"),model="within")
## Warning in pdata.frame(data, index = index, ...): duplicate couples (id-time) in resulting pdata.frame
## to find out which, use, e.g., table(index(your_pdataframe), useNA = "ifany")
summary(NE_region_model2)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = talremia_inc ~ CumPrecip + forest + mixeddev +
## TempVar + PrecipVar + AvgWinterTemp + AvgSpringPrecip + HotDryDays +
## tick_lag, data = NE_region, model = "within", index = c("county",
## "year.1"))
##
## Balanced Panel: n = 1, T = 10, N = 110
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -2.55381 -0.93210 -0.20483 0.62293 5.33474
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## CumPrecip 0.0031428 0.0012749 2.4652 0.015398 *
## forest -0.0400331 0.0075991 -5.2682 7.935e-07 ***
## mixeddev -1.8363538 0.4522750 -4.0603 9.749e-05 ***
## TempVar -0.0027142 0.0053809 -0.5044 0.615084
## PrecipVar -0.0576219 0.0173020 -3.3304 0.001216 **
## AvgWinterTemp 0.0310621 0.0498556 0.6230 0.534676
## AvgSpringPrecip -0.2959025 0.1281206 -2.3096 0.022969 *
## HotDryDays -0.0172457 0.0506623 -0.3404 0.734265
## tick_lag 0.0442506 0.0168704 2.6230 0.010080 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 419.42
## Residual Sum of Squares: 202.33
## R-Squared: 0.51759
## Adj. R-Squared: 0.47417
## F-statistic: 11.9213 on 9 and 100 DF, p-value: 1.3975e-12
#cluster robust se
library(sandwich)
clustered_se <- vcovHC(NE_region_model2, type = "HC1", cluster = "group")
coef(NE_region_model2,vcov=clustered_se)
## CumPrecip forest mixeddev TempVar PrecipVar
## 0.003142804 -0.040033130 -1.836353849 -0.002714184 -0.057621907
## AvgWinterTemp AvgSpringPrecip HotDryDays tick_lag
## 0.031062055 -0.295902452 -0.017245721 0.044250634
# Midwest region
MW_region<-project1[project1$county=="MW",]
#modelling for midwest
MW_region_model1<-plm(talremia_inc~CumPrecip+forest
+mixeddev
+TempVar
+PrecipVar
+AvgWinterTemp
+AvgSpringPrecip
+HotDryDays+tick_lag,data=MW_region,index=c("county","year"),model="within")
## Warning in pdata.frame(data, index = index, ...): duplicate couples (id-time) in resulting pdata.frame
## to find out which, use, e.g., table(index(your_pdataframe), useNA = "ifany")
summary(MW_region_model1)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = talremia_inc ~ CumPrecip + forest + mixeddev +
## TempVar + PrecipVar + AvgWinterTemp + AvgSpringPrecip + HotDryDays +
## tick_lag, data = MW_region, model = "within", index = c("county",
## "year"))
##
## Balanced Panel: n = 1, T = 10, N = 120
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -15.31231 -4.54494 -0.79275 2.99295 21.42614
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## CumPrecip -0.0325670 0.0061129 -5.3276 5.338e-07 ***
## forest 0.1275292 0.0624949 2.0406 0.043682 *
## mixeddev -10.4946630 3.7744853 -2.7804 0.006388 **
## TempVar -0.0059059 0.0161253 -0.3663 0.714882
## PrecipVar 0.4385078 0.0906899 4.8352 4.349e-06 ***
## AvgWinterTemp -0.1376964 0.2308852 -0.5964 0.552144
## AvgSpringPrecip 0.2088590 0.6059839 0.3447 0.731007
## HotDryDays 0.5832773 0.2588939 2.2530 0.026245 *
## tick_lag 0.1459472 0.0773339 1.8872 0.061765 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 10097
## Residual Sum of Squares: 5714.3
## R-Squared: 0.43405
## Adj. R-Squared: 0.38774
## F-statistic: 9.37364 on 9 and 110 DF, p-value: 1.8054e-10
#cluster robust se
library(sandwich)
clustered_se2 <- vcovHC(MW_region_model1, type = "HC1", cluster = "group")
coef(MW_region_model1,vcov=clustered_se2)
## CumPrecip forest mixeddev TempVar PrecipVar
## -0.032567047 0.127529166 -10.494662990 -0.005905907 0.438507815
## AvgWinterTemp AvgSpringPrecip HotDryDays tick_lag
## -0.137696400 0.208859043 0.583277290 0.145947197
#for southwest region
SW_region<-project1[project1$county=="SW",]
#modelling for south west region
SW_region_model1<-plm(talremia_inc~CumPrecip+forest
+mixeddev
+TempVar
+PrecipVar
+AvgWinterTemp
+AvgSpringPrecip
+HotDryDays+tick_lag,data=SW_region,index=c("county","year"),model="within")
## Warning in pdata.frame(data, index = index, ...): duplicate couples (id-time) in resulting pdata.frame
## to find out which, use, e.g., table(index(your_pdataframe), useNA = "ifany")
summary(SW_region_model1)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = talremia_inc ~ CumPrecip + forest + mixeddev +
## TempVar + PrecipVar + AvgWinterTemp + AvgSpringPrecip + HotDryDays +
## tick_lag, data = SW_region, model = "within", index = c("county",
## "year"))
##
## Balanced Panel: n = 1, T = 10, N = 40
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -14.4874 -4.7439 -2.1853 3.4764 42.0687
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## CumPrecip 0.0356011 0.0270031 1.3184 0.19735
## forest -0.0016152 0.2574322 -0.0063 0.99504
## mixeddev -7.7419448 3.9550375 -1.9575 0.05965 .
## TempVar 0.0185897 0.0678414 0.2740 0.78595
## PrecipVar -0.4546685 0.3077552 -1.4774 0.15000
## AvgWinterTemp 0.6230089 0.7871156 0.7915 0.43486
## AvgSpringPrecip 0.7133847 2.7728912 0.2573 0.79873
## HotDryDays 0.1619643 0.3603876 0.4494 0.65636
## tick_lag 0.3604576 0.3205911 1.1244 0.26978
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 6363.6
## Residual Sum of Squares: 3473.2
## R-Squared: 0.4542
## Adj. R-Squared: 0.29047
## F-statistic: 2.77396 on 9 and 30 DF, p-value: 0.017186
#cluster robust se
library(sandwich)
clustered_se3 <- vcovHC(SW_region_model1, type = "HC1", cluster = "group")
coef(SW_region_model1,vcov=clustered_se3)
## CumPrecip forest mixeddev TempVar PrecipVar
## 0.035601074 -0.001615153 -7.741944812 0.018589662 -0.454668536
## AvgWinterTemp AvgSpringPrecip HotDryDays tick_lag
## 0.623008858 0.713384724 0.161964296 0.360457623
#pacific region
pa_region<-project1[project1$county=="P",]
#modelling for pacific region
pa_region_model1<-plm(talremia_inc~CumPrecip+forest
+mixeddev
+TempVar
+PrecipVar
+AvgWinterTemp
+AvgSpringPrecip
+HotDryDays+tick_lag,data=pa_region,index=c("county","year"),model="within")
## Warning in pdata.frame(data, index = index, ...): duplicate couples (id-time) in resulting pdata.frame
## to find out which, use, e.g., table(index(your_pdataframe), useNA = "ifany")
summary(pa_region_model1)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = talremia_inc ~ CumPrecip + forest + mixeddev +
## TempVar + PrecipVar + AvgWinterTemp + AvgSpringPrecip + HotDryDays +
## tick_lag, data = pa_region, model = "within", index = c("county",
## "year"))
##
## Balanced Panel: n = 1, T = 10, N = 70
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -5.62054 -1.14415 -0.28157 1.04382 10.96404
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## CumPrecip -0.0046836 0.0023787 -1.9689 0.05358 .
## forest 0.0477614 0.0653600 0.7307 0.46778
## mixeddev -0.9603335 7.3031866 -0.1315 0.89582
## TempVar -0.0041440 0.0097092 -0.4268 0.67104
## PrecipVar 0.0588792 0.0405487 1.4521 0.15169
## AvgWinterTemp -0.0379396 0.0461967 -0.8213 0.41475
## AvgSpringPrecip 0.7882544 0.3716843 2.1208 0.03808 *
## HotDryDays 0.4970523 0.2363859 2.1027 0.03969 *
## tick_lag 0.0668377 0.0129537 5.1597 2.939e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 653.37
## Residual Sum of Squares: 320.93
## R-Squared: 0.50881
## Adj. R-Squared: 0.43513
## F-statistic: 6.90574 on 9 and 60 DF, p-value: 9.3819e-07
#cluster robust se
library(sandwich)
clustered_se4 <- vcovHC(pa_region_model1, type = "HC1", cluster = "group")
coef(pa_region_model1,vcov=clustered_se4)
## CumPrecip forest mixeddev TempVar PrecipVar
## -0.004683571 0.047761386 -0.960333546 -0.004144015 0.058879215
## AvgWinterTemp AvgSpringPrecip HotDryDays tick_lag
## -0.037939638 0.788254376 0.497052270 0.066837729
#for pacific south
ps_region<-project1[project1$county=="PS",]
#modelling for pacific south
ps_region_model1<-plm(talremia_inc~CumPrecip+forest
+mixeddev
+TempVar
+PrecipVar
+AvgWinterTemp
+AvgSpringPrecip
+HotDryDays+tick_lag
,data=ps_region,index=c("county","year"),model="within")
## Warning in pdata.frame(data, index = index, ...): duplicate couples (id-time) in resulting pdata.frame
## to find out which, use, e.g., table(index(your_pdataframe), useNA = "ifany")
summary(ps_region_model1)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = talremia_inc ~ CumPrecip + forest + mixeddev +
## TempVar + PrecipVar + AvgWinterTemp + AvgSpringPrecip + HotDryDays +
## tick_lag, data = ps_region, model = "within", index = c("county",
## "year"))
##
## Balanced Panel: n = 1, T = 10, N = 40
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -2.16519 -1.05391 -0.30397 1.02710 3.92342
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## CumPrecip 5.9123e-04 3.4276e-03 0.1725 0.86421
## forest 7.7429e-02 3.0810e-02 2.5132 0.01757 *
## mixeddev -3.7896e+01 2.0450e+01 -1.8531 0.07372 .
## TempVar -1.4731e-02 1.3943e-02 -1.0566 0.29914
## PrecipVar -5.3794e-02 3.4936e-02 -1.5398 0.13410
## AvgWinterTemp 2.3126e-02 7.6478e-02 0.3024 0.76444
## AvgSpringPrecip -6.8617e-01 4.9118e-01 -1.3970 0.17267
## HotDryDays 1.4974e-01 9.1565e-02 1.6353 0.11244
## tick_lag -2.2004e-02 2.9244e-02 -0.7524 0.45766
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 146.4
## Residual Sum of Squares: 85.756
## R-Squared: 0.41424
## Adj. R-Squared: 0.23851
## F-statistic: 2.35724 on 9 and 30 DF, p-value: 0.037781
#cluster robust se
library(sandwich)
clustered_se5 <- vcovHC(ps_region_model1, type = "HC1", cluster = "group")
coef(ps_region_model1,vcov=clustered_se5)
## CumPrecip forest mixeddev TempVar PrecipVar
## 5.912342e-04 7.742941e-02 -3.789604e+01 -1.473131e-02 -5.379431e-02
## AvgWinterTemp AvgSpringPrecip HotDryDays tick_lag
## 2.312614e-02 -6.861673e-01 1.497359e-01 -2.200407e-02
#for southeast region
SE_region<-project1[project1$county=="SE",]
#modelling for south east region
SE_region_model1<-plm(talremia_inc~CumPrecip+forest
+mixeddev
+TempVar
+PrecipVar
+AvgWinterTemp
+AvgSpringPrecip
+HotDryDays+tick_lag
,data=SE_region,index=c("county","year"),model="within")
## Warning in pdata.frame(data, index = index, ...): duplicate couples (id-time) in resulting pdata.frame
## to find out which, use, e.g., table(index(your_pdataframe), useNA = "ifany")
summary(SE_region_model1)
## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = talremia_inc ~ CumPrecip + forest + mixeddev +
## TempVar + PrecipVar + AvgWinterTemp + AvgSpringPrecip + HotDryDays +
## tick_lag, data = SE_region, model = "within", index = c("county",
## "year"))
##
## Balanced Panel: n = 1, T = 10, N = 120
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -20.80932 -4.20306 -0.55705 2.63197 74.28454
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## CumPrecip -0.0022348 0.0094197 -0.2373 0.8129
## forest -0.2752164 0.0573016 -4.8029 4.971e-06 ***
## mixeddev 0.3240405 0.5238704 0.6186 0.5375
## TempVar 0.0554697 0.0383021 1.4482 0.1504
## PrecipVar -0.0708119 0.1177459 -0.6014 0.5488
## AvgWinterTemp 0.0020237 0.3230981 0.0063 0.9950
## AvgSpringPrecip 1.0073852 0.8419796 1.1964 0.2341
## HotDryDays 0.0655043 0.1967337 0.3330 0.7398
## tick_lag 0.0693735 0.1104660 0.6280 0.5313
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 18682
## Residual Sum of Squares: 12813
## R-Squared: 0.31414
## Adj. R-Squared: 0.25802
## F-statistic: 5.59807 on 9 and 110 DF, p-value: 2.399e-06
#cluster robust se
library(sandwich)
clustered_se6 <- vcovHC(SE_region_model1, type = "HC1", cluster = "group")
coef(SE_region_model1,vcov=clustered_se6)
## CumPrecip forest mixeddev TempVar PrecipVar
## -0.002234845 -0.275216448 0.324040509 0.055469680 -0.070811905
## AvgWinterTemp AvgSpringPrecip HotDryDays tick_lag
## 0.002023730 1.007385214 0.065504280 0.069373538