R Notebook for Lab 10 Space time models
Reading in the data
library(sf) ## Warning: package 'sf' was built under R version 4.4.2## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE## Data management
setwd("C:/Users/ssrini06/Box/SpatStat_Labs/Assignment10_SpTime/Lab10_Montana")
fire_notclean <- read_sf("Montana_fires_proj_corr.shp")
# Data frame for the table with subset of the variables remove NA
library("tidyr")
fire_clean <- na.omit(fire_notclean[c(2:4,11:18, 20:29,37:65)])
#Make sure date is read as a date
class(fire_clean$FireDiscov)## [1] "Date"#if you need to convert then use this code
fire_clean$NewDate <- as.Date(fire_clean$FireDiscov,"%m-%d-%y")
fire <- fire_clean[order(as.Date(fire_clean$NewDate, format="%m/%d/%Y")),]
fire_table <- data.frame(fire)
#removing NA locations and getting only the X Y coordinate and date
fire_xyt <- na.omit(fire_table[c(49:50,6,1)])
fire_xy <- na.omit(fire_table[c(49:50)])Exploratory space time data analysis
#Create a new variable YearCat for time intervals: five categories
fire$YearCat4 <- "1"
fire$YearCat4[fire_table$NewDate > as.Date("2021-01-01") & fire_table$NewDate < as.Date("2022-01-01")] <- "2"
fire$YearCat4[fire_table$NewDate > as.Date("2022-01-01") & fire_table$NewDate < as.Date("2023-01-01")] <- "3"
fire$YearCat4[fire_table$NewDate > as.Date("2023-01-01") & fire_table$NewDate < as.Date("2024-01-01")] <- "4"
summary.factor(fire$YearCat4)## 1 2 3 4
## 193 239 178 132#two categories: before and after Jan 1 22
fire$YearCat <- "1"
fire$YearCat[fire_table$NewDate > as.Date("2022-01-01")] <- "2"
summary.factor(fire$YearCat)## 1 2
## 432 310#https://data-nifc.opendata.arcgis.com/datasets/nifc::inform-fire-occurrence-data-records/about
summary.factor(fire$Predominan)## Brush Grass Grass-Shrub Nonburnable Slash Timber
## 22 140 92 3 39 446summary.factor(fire$POOCounty)## Beaverhead Big Horn Blaine Broadwater Carbon
## 23 12 3 5 7
## Carter Cascade Custer Dawson Deer Lodge
## 22 5 11 2 6
## Fergus Flathead Gallatin Garfield Glacier
## 4 63 20 11 1
## Granite Jefferson Judith Basin Lake Lewis and Clark
## 12 25 4 24 29
## Liberty Lincoln Madison McCone Meagher
## 1 8 4 8 5
## Mineral Missoula Musselshell Park Petroleum
## 36 94 1 21 3
## Phillips Pondera Powder River Powell Prairie
## 2 1 81 19 6
## Ravalli Roosevelt Rosebud Sanders Silver Bow
## 54 2 26 44 20
## Stillwater Sweet Grass Teton Wheatland Wibaux
## 2 6 3 1 3
## Yellowstone
## 2#Days that fire was not contained
summary(fire$FireDays)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 1.00 8.00 19.75 27.00 155.00#Acres of fire
summary(fire$Calculated)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.100 0.250 401.263 6.968 29884.810# Exploratory Temporal Data Analysis
library(ggplot2)
#changing from e notation
options(scipen=999)
p <- ggplot(fire_table, aes(x = NewDate, y = Calculated))
p + geom_line() + geom_point() + xlab("Years") + ylab("Acreage") 
p + geom_line(aes(color = Predominan)) + geom_point() + xlab("Time") + ylab("Acreage") 
# Plot the points
# change the basemap to ESRI Topomap
plot(st_geometry(fire))
library(RColorBrewer)
library(tmap)## Warning: package 'tmap' was built under R version 4.4.3tm_shape(fire) +
tm_dots(fill="Calculated",
fill.scale = tm_scale(values="brewer.reds", style = "quantile"),
fill.legend = tm_legend(title = "Acreage", size = 0.8))+
tmap_mode("view") ## ℹ tmap mode set to "view".
## Registered S3 method overwritten by 'jsonify':
## method from
## print.json jsonlite#plot the time line
library(spacetime)## Warning: package 'spacetime' was built under R version 4.4.3library(sp)
crs_fire <- proj4string(as(fire, 'Spatial'))
pts <- SpatialPoints(fire_xy, CRS("+init=epsg:32100"))## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.day <- as.Date(fire_xyt$FireDiscov)
stiObj <- STI(pts, day)
plot(stiObj)
# plot shows the space-time (time=x, space=y)
#a lot of points at a location in the plot
#suggests Space and Time clustering
class(stiObj)## [1] "STI"
## attr(,"package")
## [1] "spacetime"#STIDF is a space time data frame
spat_part <- SpatialPoints(fire_xyt[,1:2])
temp_part <- fire_xyt[,3]
data_part <- fire_xyt[,4]
STIDFObj <- STIDF(sp = spat_part,
time = temp_part,
data = data.frame(data_part))
#plotting the space time data
library(RColorBrewer)
stplot(STIDFObj, number=8, col.regions = brewer.pal(7, "YlOrRd"),cuts=7)
Section 2.1 Space time descriptive statistics Knox and Mantel tests
# install.packages("surveillance")
library(surveillance)## Warning: package 'surveillance' was built under R version 4.4.2## Loading required package: xtable## This is surveillance 1.24.1; see 'package?surveillance' or
## https://surveillance.R-Forge.R-project.org/ for an overview.## Perform the Knox test using the Poisson approximation
#for this data:
#dt is time step in days
#ds is distance in space in meters
knoxtest <- knox(
dt = dist(fire$FireDiscov), eps.t = 15,
ds = dist(st_coordinates(fire)), eps.s = 50000,
simulate.p.value = FALSE
)
knoxtest##
## Knox test with Poisson approximation
##
## data: dt = dist(fire$FireDiscov) and ds = dist(st_coordinates(fire))
## number of close pairs = 1473, lambda = 1091.4, p-value <
## 0.00000000000000022
## alternative hypothesis: true number is greater than 1091.419
##
## contingency table:
## ds
## dt <= 50000 > 50000
## <= 15 1473 16160
## > 15 15543 241735## Obtain the p-value via a Monte Carlo permutation test
# B is the number of permutations
# Higher number will take more time
knoxtest <- knox(
dt = dist(fire$FireDiscov), eps.t = 14,
ds = dist(st_coordinates(fire)), eps.s = 50000,
simulate.p.value = TRUE, B=100,
.parallel = 2, .seed = 1, .verbose = FALSE
)
#contingency table compares all point pairs (n*n-1)/2
knoxtest##
## Knox test with simulated p-value
##
## data: dt = dist(fire$FireDiscov) and ds = dist(st_coordinates(fire))
## number of close pairs = 1420, B = 100, p-value = 0.009901
## alternative hypothesis: true number is greater than 1033.732
##
## contingency table:
## ds
## dt <= 50000 > 50000
## <= 14 1420 15281
## > 14 15596 242614#install.packages("ade4")
library(ade4)## Warning: package 'ade4' was built under R version 4.4.3##
## Attaching package: 'ade4'## The following object is masked from 'package:surveillance':
##
## scoredt <- dist(fire_xyt$FireDiscov)
#distances are how far apart they are in time days
# five categories
as.matrix(dt)[1:5, 1:5]## 1 2 3 4 5
## 1 0 58 64 72 77
## 2 58 0 6 14 19
## 3 64 6 0 8 13
## 4 72 14 8 0 5
## 5 77 19 13 5 0ds <- dist(coordinates(SpatialPoints(fire_xyt[,1:2])))
#distances are actual distances in meters
#five categories
as.matrix(ds)[1:5, 1:5]## 1 2 3 4 5
## 1 0.00 73123.44 31060.90 62842.34 66711.83
## 2 73123.44 0.00 97051.33 127499.51 6508.06
## 3 31060.90 97051.33 0.00 32123.80 90560.88
## 4 62842.34 127499.51 32123.80 0.00 121082.50
## 5 66711.83 6508.06 90560.88 121082.50 0.00mantel.rtest(ds, dt, nrepet = 499)## Warning in is.euclid(m2): Zero distance(s)## Monte-Carlo test
## Call: mantel.rtest(m1 = ds, m2 = dt, nrepet = 499)
##
## Observation: 0.01778736
##
## Based on 499 replicates
## Simulated p-value: 0.03
## Alternative hypothesis: greater
##
## Std.Obs Expectation Variance
## 1.86641262205 0.00044336004 0.00008635415Section 2.2 Space time K cross tests
# point patterns
library(spatstat)## Warning: package 'spatstat' was built under R version 4.4.3## Loading required package: spatstat.data## Warning: package 'spatstat.data' was built under R version 4.4.3## Loading required package: spatstat.univar## Warning: package 'spatstat.univar' was built under R version 4.4.3## spatstat.univar 3.1-2## Loading required package: spatstat.geom## Warning: package 'spatstat.geom' was built under R version 4.4.3## spatstat.geom 3.3-6##
## Attaching package: 'spatstat.geom'## The following object is masked from 'package:ade4':
##
## disc## Loading required package: spatstat.random## Warning: package 'spatstat.random' was built under R version 4.4.3## spatstat.random 3.3-2## Loading required package: spatstat.explore## Warning: package 'spatstat.explore' was built under R version 4.4.2## Loading required package: nlme## spatstat.explore 3.3-4## Loading required package: spatstat.model## Warning: package 'spatstat.model' was built under R version 4.4.2## Loading required package: rpart## spatstat.model 3.3-4## Loading required package: spatstat.linnet## Warning: package 'spatstat.linnet' was built under R version 4.4.3## spatstat.linnet 3.2-5##
## spatstat 3.3-1
## For an introduction to spatstat, type 'beginner'fire_sp <- as_Spatial(fire)
firepp <- as.ppp(fire)
firepp.km <- rescale(firepp, 1000, "km")
window <- as.owin(firepp.km)
#marks are the types of points being compared
# pre 2022 and post
marks(firepp.km) <- factor(fire$YearCat)
#could use inhomogeneous because the points are
# not homogeneously distributed
#comparing point pattern of two different years
#Kcross suggests points in both time periods cluster
# with each other until about 50km
kcrossresults <- Kcross(firepp.km)
plot(kcrossresults)
#Compare multiple years
#if its above the poisson more clustered than expected
#along the diagonal its a K not a K cross
marks(firepp.km) <- factor(fire$YearCat4)
plot(alltypes(firepp.km, Kcross, envelope=TRUE))
#for every time period the fires cluster
#they are close to fires from other time periodsSection 3.0 Kriging in space and time
#empirical variogram
#takes some time
library(gstat)## Warning: package 'gstat' was built under R version 4.4.3##
## Attaching package: 'gstat'## The following object is masked from 'package:spatstat.explore':
##
## idwvar <- variogramST(data_part~1,data=STIDFObj, tlags=1:10,tunit="weeks",assumeRegular=F,na.omit=T)
# interpret like usual variogram which lag works best?
#For higher distances the differences in Area should be higher
# a lag of 1 week seems closest to an expected variogram
plot(var, map=F)
#for lag of 1-3 weeks you see increasing semivariance with distance
plot(var, map=T)
plot(var, wireframe=T)
#demo shows you how to space time krige for a dataset for pollution # where there is a space time trend unlike our data for fires
#demo(stkrige) 4.1 Regressions with time lags
#removing NA locations
fire_model_data <- data.frame(na.omit(fire))
library(stargazer)##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazer#save summary to text or html file
stargazer(fire_model_data,type = "text", out="summary.txt")##
## ================================================================
## Statistic N Mean St. Dev. Min Max
## ----------------------------------------------------------------
## Calculated 742 401.263 2,296.127 0.000 29,884.810
## FinalAcres 742 0.000 0.000 0 0
## InitialLon 742 -105.164 25.684 -125.700 0.000
## IsFireCaus 742 0.152 0.360 0 1
## IsFSAssist 742 0.799 0.401 0 1
## IsReimburs 742 0.000 0.000 0 0
## IsTrespass 742 0.062 0.241 0 1
## PercentCon 742 0.000 0.000 0 0
## FireCauseI 742 2.049 0.965 1 3
## FireCauseG 742 4.491 4.372 1 13
## FireCauseS 742 7.987 18.621 0 77
## FireCaus_1 742 0.488 4.505 0 60
## FireCauseA 742 0.071 0.489 0 4
## FireCaus_2 742 0.515 2.896 0 19
## IsCausePro 742 0.011 0.103 0 1
## AcresBIA 742 6.116 118.223 0.000 2,939.440
## AcresBLM 742 40.943 509.365 0.000 10,870.600
## AcresBOR 742 0.000 0.000 0 0
## AcresDOD 742 3.051 80.736 0.000 2,198.390
## AcresDOE 742 0.000 0.000 0 0
## AcresNPS 742 2.585 69.108 0.000 1,882.310
## AcresUSFS 742 191.794 1,458.324 0.000 23,953.170
## AcresUSFWS 742 8.277 214.758 0.000 5,847.680
## AcresForei 742 0.000 0.000 0 0
## AcresTriba 742 75.490 939.879 0.000 18,806.970
## AcresCity 742 0.000 0.000 0 0
## AcresCount 742 0.000 0.000 0 0
## AcresState 742 8.740 95.870 0.000 1,566.200
## AcresPriva 742 50.892 502.707 0.000 10,273.070
## AcresANCSA 742 0.000 0.000 0 0
## AcresOther 742 0.000 0.000 0 0
## AcresOth_1 742 0.001 0.018 0.000 0.500
## FireDays 742 19.747 28.275 0 155
## XCoord 742 453,684.700 266,663.400 122,933.800 1,020,310.000
## YCoord 742 268,805.000 113,843.500 41,429.320 534,360.900
## ----------------------------------------------------------------stargazer(fire_model_data,type = "html", out="summary.html")##
## <table style="text-align:center"><tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Statistic</td><td>N</td><td>Mean</td><td>St. Dev.</td><td>Min</td><td>Max</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Calculated</td><td>742</td><td>401.263</td><td>2,296.127</td><td>0.000</td><td>29,884.810</td></tr>
## <tr><td style="text-align:left">FinalAcres</td><td>742</td><td>0.000</td><td>0.000</td><td>0</td><td>0</td></tr>
## <tr><td style="text-align:left">InitialLon</td><td>742</td><td>-105.164</td><td>25.684</td><td>-125.700</td><td>0.000</td></tr>
## <tr><td style="text-align:left">IsFireCaus</td><td>742</td><td>0.152</td><td>0.360</td><td>0</td><td>1</td></tr>
## <tr><td style="text-align:left">IsFSAssist</td><td>742</td><td>0.799</td><td>0.401</td><td>0</td><td>1</td></tr>
## <tr><td style="text-align:left">IsReimburs</td><td>742</td><td>0.000</td><td>0.000</td><td>0</td><td>0</td></tr>
## <tr><td style="text-align:left">IsTrespass</td><td>742</td><td>0.062</td><td>0.241</td><td>0</td><td>1</td></tr>
## <tr><td style="text-align:left">PercentCon</td><td>742</td><td>0.000</td><td>0.000</td><td>0</td><td>0</td></tr>
## <tr><td style="text-align:left">FireCauseI</td><td>742</td><td>2.049</td><td>0.965</td><td>1</td><td>3</td></tr>
## <tr><td style="text-align:left">FireCauseG</td><td>742</td><td>4.491</td><td>4.372</td><td>1</td><td>13</td></tr>
## <tr><td style="text-align:left">FireCauseS</td><td>742</td><td>7.987</td><td>18.621</td><td>0</td><td>77</td></tr>
## <tr><td style="text-align:left">FireCaus_1</td><td>742</td><td>0.488</td><td>4.505</td><td>0</td><td>60</td></tr>
## <tr><td style="text-align:left">FireCauseA</td><td>742</td><td>0.071</td><td>0.489</td><td>0</td><td>4</td></tr>
## <tr><td style="text-align:left">FireCaus_2</td><td>742</td><td>0.515</td><td>2.896</td><td>0</td><td>19</td></tr>
## <tr><td style="text-align:left">IsCausePro</td><td>742</td><td>0.011</td><td>0.103</td><td>0</td><td>1</td></tr>
## <tr><td style="text-align:left">AcresBIA</td><td>742</td><td>6.116</td><td>118.223</td><td>0.000</td><td>2,939.440</td></tr>
## <tr><td style="text-align:left">AcresBLM</td><td>742</td><td>40.943</td><td>509.365</td><td>0.000</td><td>10,870.600</td></tr>
## <tr><td style="text-align:left">AcresBOR</td><td>742</td><td>0.000</td><td>0.000</td><td>0</td><td>0</td></tr>
## <tr><td style="text-align:left">AcresDOD</td><td>742</td><td>3.051</td><td>80.736</td><td>0.000</td><td>2,198.390</td></tr>
## <tr><td style="text-align:left">AcresDOE</td><td>742</td><td>0.000</td><td>0.000</td><td>0</td><td>0</td></tr>
## <tr><td style="text-align:left">AcresNPS</td><td>742</td><td>2.585</td><td>69.108</td><td>0.000</td><td>1,882.310</td></tr>
## <tr><td style="text-align:left">AcresUSFS</td><td>742</td><td>191.794</td><td>1,458.324</td><td>0.000</td><td>23,953.170</td></tr>
## <tr><td style="text-align:left">AcresUSFWS</td><td>742</td><td>8.277</td><td>214.758</td><td>0.000</td><td>5,847.680</td></tr>
## <tr><td style="text-align:left">AcresForei</td><td>742</td><td>0.000</td><td>0.000</td><td>0</td><td>0</td></tr>
## <tr><td style="text-align:left">AcresTriba</td><td>742</td><td>75.490</td><td>939.879</td><td>0.000</td><td>18,806.970</td></tr>
## <tr><td style="text-align:left">AcresCity</td><td>742</td><td>0.000</td><td>0.000</td><td>0</td><td>0</td></tr>
## <tr><td style="text-align:left">AcresCount</td><td>742</td><td>0.000</td><td>0.000</td><td>0</td><td>0</td></tr>
## <tr><td style="text-align:left">AcresState</td><td>742</td><td>8.740</td><td>95.870</td><td>0.000</td><td>1,566.200</td></tr>
## <tr><td style="text-align:left">AcresPriva</td><td>742</td><td>50.892</td><td>502.707</td><td>0.000</td><td>10,273.070</td></tr>
## <tr><td style="text-align:left">AcresANCSA</td><td>742</td><td>0.000</td><td>0.000</td><td>0</td><td>0</td></tr>
## <tr><td style="text-align:left">AcresOther</td><td>742</td><td>0.000</td><td>0.000</td><td>0</td><td>0</td></tr>
## <tr><td style="text-align:left">AcresOth_1</td><td>742</td><td>0.001</td><td>0.018</td><td>0.000</td><td>0.500</td></tr>
## <tr><td style="text-align:left">FireDays</td><td>742</td><td>19.747</td><td>28.275</td><td>0</td><td>155</td></tr>
## <tr><td style="text-align:left">XCoord</td><td>742</td><td>453,684.700</td><td>266,663.400</td><td>122,933.800</td><td>1,020,310.000</td></tr>
## <tr><td style="text-align:left">YCoord</td><td>742</td><td>268,805.000</td><td>113,843.500</td><td>41,429.320</td><td>534,360.900</td></tr>
## <tr><td colspan="6" style="border-bottom: 1px solid black"></td></tr></table>#Model with all the data
#outcome is FireDays and explanatory variable is Area
ols <- lm(FireDays ~ Calculated, data=fire_model_data)
summary(ols)##
## Call:
## lm(formula = FireDays ~ Calculated, data = fire_model_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -87.51 -16.67 -10.67 8.33 135.33
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.6695435 0.9567848 18.47 <0.0000000000000002 ***
## Calculated 0.0051764 0.0004107 12.60 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 25.67 on 740 degrees of freedom
## Multiple R-squared: 0.1767, Adjusted R-squared: 0.1756
## F-statistic: 158.8 on 1 and 740 DF, p-value: < 0.00000000000000022# Using time series data with quarters fire_ts
library(zoo)## Warning: package 'zoo' was built under R version 4.4.3##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericfire_model_data$fire_qts <- as.yearqtr(fire_model_data$NewDate, format = "%Y-%m-%d")
head(fire_model_data) # Aggregate data to quarters
# Making sure it reads dplyr function summarize and not from Hmisc
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:nlme':
##
## collapse## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionfire_ts <- fire_model_data %>% group_by(fire_qts) %>%
dplyr::summarize(Days_SUM = sum(FireDays), ACRE_SUM= sum(Calculated))
head(fire_ts) # MOdel with quarters instead of all points
ols2 <- lm(Days_SUM ~ ACRE_SUM, data=fire_ts)
summary(ols2)##
## Call:
## lm(formula = Days_SUM ~ ACRE_SUM, data = fire_ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1164.70 -78.00 -21.14 154.58 1011.93
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 133.615152 158.818595 0.841 0.415
## ACRE_SUM 0.042480 0.004412 9.628 0.000000279 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 513.1 on 13 degrees of freedom
## Multiple R-squared: 0.877, Adjusted R-squared: 0.8676
## F-statistic: 92.71 on 1 and 13 DF, p-value: 0.0000002787library(Hmisc)## Warning: package 'Hmisc' was built under R version 4.4.3##
## Attaching package: 'Hmisc'## The following objects are masked from 'package:dplyr':
##
## src, summarize## The following objects are masked from 'package:xtable':
##
## label, label<-## The following objects are masked from 'package:base':
##
## format.pval, units#using a lag for area 1 day or 7 day
fire_model_data$lag1time <- Lag(fire_model_data$FireDays,1)
fire_model_data$lag7time <- Lag(fire_model_data$FireDays,7)
olstime <- lm(FireDays ~ Calculated + Lag(FireDays,1), data=fire_model_data)
summary(olstime)##
## Call:
## lm(formula = FireDays ~ Calculated + Lag(FireDays, 1), data = fire_model_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -95.150 -15.845 -9.847 7.809 137.154
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.8445255 1.1539632 13.731 < 0.0000000000000002 ***
## Calculated 0.0051270 0.0004093 12.526 < 0.0000000000000002 ***
## Lag(FireDays, 1) 0.0944767 0.0332428 2.842 0.00461 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 25.56 on 738 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1856, Adjusted R-squared: 0.1834
## F-statistic: 84.09 on 2 and 738 DF, p-value: < 0.00000000000000022#quarter based lag
olstime2 <- lm(Days_SUM ~ ACRE_SUM + Lag(Days_SUM,1), data=fire_ts)
summary(olstime2)##
## Call:
## lm(formula = Days_SUM ~ ACRE_SUM + Lag(Days_SUM, 1), data = fire_ts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1164.59 -93.54 1.44 184.76 1022.00
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 109.387141 229.691023 0.476 0.643
## ACRE_SUM 0.042705 0.005013 8.520 0.00000357 ***
## Lag(Days_SUM, 1) 0.016775 0.110646 0.152 0.882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 557.2 on 11 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.8741, Adjusted R-squared: 0.8512
## F-statistic: 38.18 on 2 and 11 DF, p-value: 0.00001124#saving results table
stargazer(ols,olstime, ols2, olstime2, column.labels = c("All data", "All data", "Quarters Data", "Quarters Data"), digits=1, type = "html", out="ols.html")##
## <table style="text-align:center"><tr><td colspan="5" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="4"><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="4" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td colspan="2">FireDays</td><td colspan="2">Days_SUM</td></tr>
## <tr><td style="text-align:left"></td><td>All data</td><td>All data</td><td>Quarters Data</td><td>Quarters Data</td></tr>
## <tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td><td>(3)</td><td>(4)</td></tr>
## <tr><td colspan="5" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Calculated</td><td>0.01<sup>***</sup></td><td>0.01<sup>***</sup></td><td></td><td></td></tr>
## <tr><td style="text-align:left"></td><td>(0.000)</td><td>(0.000)</td><td></td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Lag(FireDays, 1)</td><td></td><td>0.1<sup>***</sup></td><td></td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td>(0.03)</td><td></td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">ACRE_SUM</td><td></td><td></td><td>0.04<sup>***</sup></td><td>0.04<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td>(0.004)</td><td>(0.01)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Lag(Days_SUM, 1)</td><td></td><td></td><td></td><td>0.02</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td>(0.1)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>17.7<sup>***</sup></td><td>15.8<sup>***</sup></td><td>133.6</td><td>109.4</td></tr>
## <tr><td style="text-align:left"></td><td>(1.0)</td><td>(1.2)</td><td>(158.8)</td><td>(229.7)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td></tr>
## <tr><td colspan="5" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>742</td><td>741</td><td>15</td><td>14</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.2</td><td>0.2</td><td>0.9</td><td>0.9</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.2</td><td>0.2</td><td>0.9</td><td>0.9</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>25.7 (df = 740)</td><td>25.6 (df = 738)</td><td>513.1 (df = 13)</td><td>557.2 (df = 11)</td></tr>
## <tr><td style="text-align:left">F Statistic</td><td>158.8<sup>***</sup> (df = 1; 740)</td><td>84.1<sup>***</sup> (df = 2; 738)</td><td>92.7<sup>***</sup> (df = 1; 13)</td><td>38.2<sup>***</sup> (df = 2; 11)</td></tr>
## <tr><td colspan="5" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="4" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>4.2 Regressions with with fixed and random effects for time
#panel regression with fixed effects for year
require(plm)## Loading required package: plm## Warning: package 'plm' was built under R version 4.4.2##
## Attaching package: 'plm'## The following objects are masked from 'package:dplyr':
##
## between, lag, leadfire_model_data$FIRE_YEAR <- format(fire_model_data$NewDate, format="%Y")
fixed_effects_fire <-plm(FireDays ~ Calculated,
index = "FIRE_YEAR",
model="within",
data=fire_model_data)
summary(fixed_effects_fire)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = FireDays ~ Calculated, data = fire_model_data,
## model = "within", index = "FIRE_YEAR")
##
## Unbalanced Panel: n = 4, T = 132-239, N = 742
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -90.5645 -15.1258 -9.1253 6.6104 129.9676
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## Calculated 0.00505327 0.00040478 12.484 < 0.00000000000000022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 568460
## Residual Sum of Squares: 469230
## R-Squared: 0.17455
## Adj. R-Squared: 0.17007
## F-statistic: 155.849 on 1 and 737 DF, p-value: < 0.000000000000000222# n = number of groups/panels, T = number of time periods, N = total number of observations
# with dummies for time intervals
fixed_effects_fireqt <-plm(FireDays ~ Calculated,
index = "YearCat4",
model="within",
data=fire_model_data)
summary(fixed_effects_fireqt)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = FireDays ~ Calculated, data = fire_model_data,
## model = "within", index = "YearCat4")
##
## Unbalanced Panel: n = 4, T = 132-239, N = 742
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -90.5645 -15.1258 -9.1253 6.6104 129.9676
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## Calculated 0.00505327 0.00040478 12.484 < 0.00000000000000022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 568460
## Residual Sum of Squares: 469230
## R-Squared: 0.17455
## Adj. R-Squared: 0.17007
## F-statistic: 155.849 on 1 and 737 DF, p-value: < 0.000000000000000222# Extracting fixed effects with fixef
# shows how many days of fire increases (or not) for unit acreage increase
# in each time period
# note that only some years are significant and can be negative
fixef(fixed_effects_fire, type = "dmean")## 2020 2021 2022 2023
## -7.4400 5.3122 2.1165 -1.5941summary(fixef(fixed_effects_fire, type = "dmean"))## Estimate Std. Error t-value Pr(>|t|)
## 2020 -7.4400 1.8177 -4.0930 0.00004728 ***
## 2021 5.3122 1.6489 3.2216 0.001331 **
## 2022 2.1165 1.8953 1.1167 0.264475
## 2023 -1.5941 2.2067 -0.7224 0.470283
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fixef(fixed_effects_fireqt, type = "dmean")## 1 2 3 4
## -7.4400 5.3122 2.1165 -1.5941summary(fixef(fixed_effects_fireqt, type = "dmean"))## Estimate Std. Error t-value Pr(>|t|)
## 1 -7.4400 1.8177 -4.0930 0.00004728 ***
## 2 5.3122 1.6489 3.2216 0.001331 **
## 3 2.1165 1.8953 1.1167 0.264475
## 4 -1.5941 2.2067 -0.7224 0.470283
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# if p is small used fixed effects
# instead of ols
pFtest(fixed_effects_fire, ols)##
## F test for individual effects
##
## data: FireDays ~ Calculated
## F = 9.6864, df1 = 3, df2 = 737, p-value = 0.000002824
## alternative hypothesis: significant effectspFtest(fixed_effects_fireqt, ols)##
## F test for individual effects
##
## data: FireDays ~ Calculated
## F = 9.6864, df1 = 3, df2 = 737, p-value = 0.000002824
## alternative hypothesis: significant effects#random effects
random_effects_fire <-plm(FireDays ~ Calculated,
index = c("FIRE_YEAR"),
model="random",
data=fire_model_data)
summary(random_effects_fire)## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = FireDays ~ Calculated, data = fire_model_data,
## model = "random", index = c("FIRE_YEAR"))
##
## Unbalanced Panel: n = 4, T = 132-239, N = 742
##
## Effects:
## var std.dev share
## idiosyncratic 636.678 25.232 0.973
## individual 17.947 4.236 0.027
## theta:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.5398 0.5924 0.6060 0.6020 0.6405 0.6405
##
## Residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -88.843 -14.726 -10.071 0.066 7.173 132.010
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) 17.36726710 2.32654593 7.4648 0.00000000000008341 ***
## Calculated 0.00507043 0.00040496 12.5208 < 0.00000000000000022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 571490
## Residual Sum of Squares: 471950
## R-Squared: 0.17418
## Adj. R-Squared: 0.17306
## Chisq: 156.769 on 1 DF, p-value: < 0.000000000000000222random.effects(random_effects_fire, type = "dmean")## 2020 2021 2022 2023
## -5.9903837 4.9231601 2.0536784 -0.9864548random_effects_fireqt <-plm(FireDays ~ Calculated,
index = c("YearCat4"),
model="random",
data=fire_model_data)
summary(random_effects_fireqt)## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = FireDays ~ Calculated, data = fire_model_data,
## model = "random", index = c("YearCat4"))
##
## Unbalanced Panel: n = 4, T = 132-239, N = 742
##
## Effects:
## var std.dev share
## idiosyncratic 636.678 25.232 0.973
## individual 17.947 4.236 0.027
## theta:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.5398 0.5924 0.6060 0.6020 0.6405 0.6405
##
## Residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -88.843 -14.726 -10.071 0.066 7.173 132.010
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) 17.36726710 2.32654593 7.4648 0.00000000000008341 ***
## Calculated 0.00507043 0.00040496 12.5208 < 0.00000000000000022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 571490
## Residual Sum of Squares: 471950
## R-Squared: 0.17418
## Adj. R-Squared: 0.17306
## Chisq: 156.769 on 1 DF, p-value: < 0.000000000000000222random.effects(random_effects_fireqt, type = "dmean")## 1 2 3 4
## -5.9903837 4.9231601 2.0536784 -0.9864548summary(random.effects(random_effects_fireqt, type = "dmean"))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -5.9904 -2.2374 0.5336 0.0000 2.7710 4.9232phtest(fixed_effects_fire, random_effects_fire)##
## Hausman Test
##
## data: FireDays ~ Calculated
## chisq = 2.0102, df = 1, p-value = 0.1562
## alternative hypothesis: one model is inconsistentphtest(fixed_effects_fireqt, random_effects_fireqt)##
## Hausman Test
##
## data: FireDays ~ Calculated
## chisq = 2.0102, df = 1, p-value = 0.1562
## alternative hypothesis: one model is inconsistent#indicates that the null hypothesis
#that the individual random effects are exogenous is rejected
#which makes the random effects equation inconsistent
#the fixed effects model is better
stargazer(ols, fixed_effects_fire, fixed_effects_fireqt, digits=1, type = "html", out="ols_ml.html") ##
## <table style="text-align:center"><tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="3"><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="3" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td colspan="3">FireDays</td></tr>
## <tr><td style="text-align:left"></td><td><em>OLS</em></td><td colspan="2"><em>panel</em></td></tr>
## <tr><td style="text-align:left"></td><td><em></em></td><td colspan="2"><em>linear</em></td></tr>
## <tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td><td>(3)</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Calculated</td><td>0.01<sup>***</sup></td><td>0.01<sup>***</sup></td><td>0.01<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.000)</td><td>(0.000)</td><td>(0.000)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>17.7<sup>***</sup></td><td></td><td></td></tr>
## <tr><td style="text-align:left"></td><td>(1.0)</td><td></td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>742</td><td>742</td><td>742</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.2</td><td>0.2</td><td>0.2</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.2</td><td>0.2</td><td>0.2</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>25.7 (df = 740)</td><td></td><td></td></tr>
## <tr><td style="text-align:left">F Statistic</td><td>158.8<sup>***</sup> (df = 1; 740)</td><td>155.8<sup>***</sup> (df = 1; 737)</td><td>155.8<sup>***</sup> (df = 1; 737)</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="3" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>Section 4.3 Spatial regressions with time lags
library(spatialreg)## Warning: package 'spatialreg' was built under R version 4.4.3## Loading required package: spData## Warning: package 'spData' was built under R version 4.4.3## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(spdep)## Warning: package 'spdep' was built under R version 4.4.3##
## Attaching package: 'spdep'## The following objects are masked from 'package:spatialreg':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOption## The following object is masked from 'package:ade4':
##
## mstree#distance based neighbors
#generate distances between centroids for the nearest neighbors
fire_coords <- st_coordinates(fire)
fireIDs <- row.names(as(fire_sp, "data.frame"))
fire_nn1 <- knn2nb(knearneigh(fire_coords, k=1), row.names= fireIDs)## Warning in knn2nb(knearneigh(fire_coords, k = 1), row.names = fireIDs):
## neighbour object has 205 sub-graphsfiredist <-unlist(nbdists(fire_nn1, fire_coords))
#shows the min mean and max distance between centroids
summary(firedist)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.61 1848.96 4183.28 6859.38 8762.07 151829.00#distance threshold to ensure all fires have a neighbor 78km
dist_threshold <-max(firedist)
#generating distance based neighbors
fire_dist_threshold <-dnearneigh(fire_coords, d1=0,
d2=1*dist_threshold,
row.names= fireIDs)
W_dist_rowstd <- nb2listw(fire_dist_threshold, zero.policy=TRUE)
#Test to see if you need a spatial regression
lm.RStests(olstime, W_dist_rowstd, test="all", zero.policy=T)##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = FireDays ~ Calculated + Lag(FireDays, 1), data =
## fire_model_data)
## test weights: W_dist_rowstd
##
## RSerr = 68.542, df = 1, p-value < 0.00000000000000022
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = FireDays ~ Calculated + Lag(FireDays, 1), data =
## fire_model_data)
## test weights: W_dist_rowstd
##
## RSlag = 92.952, df = 1, p-value < 0.00000000000000022
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = FireDays ~ Calculated + Lag(FireDays, 1), data =
## fire_model_data)
## test weights: W_dist_rowstd
##
## adjRSerr = 2.2086, df = 1, p-value = 0.1372
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = FireDays ~ Calculated + Lag(FireDays, 1), data =
## fire_model_data)
## test weights: W_dist_rowstd
##
## adjRSlag = 26.618, df = 1, p-value = 0.0000002479
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = FireDays ~ Calculated + Lag(FireDays, 1), data =
## fire_model_data)
## test weights: W_dist_rowstd
##
## SARMA = 95.16, df = 2, p-value < 0.00000000000000022# Error regression
# If it does not run in R try it in Geoda
fire_spatial_error <- errorsarlm(FireDays ~ Calculated + Lag(FireDays,1), data=fire_model_data, W_dist_rowstd, zero.policy=T, tol.solve=20.0e-20)
summary(fire_spatial_error)##
## Call:errorsarlm(formula = FireDays ~ Calculated + Lag(FireDays, 1),
## data = fire_model_data, listw = W_dist_rowstd, zero.policy = T,
## tol.solve = 0.0000000000000000002)
##
## Residuals:
## Min 1Q Median 3Q Max
## -93.7333 -13.9446 -9.8661 7.0699 130.9382
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 17.51065268 3.52437333 4.9684 0.0000006749
## Calculated 0.00483111 0.00040013 12.0739 < 0.00000000000000022
## Lag(FireDays, 1) 0.08782480 0.03255621 2.6976 0.006983
##
## Lambda: 0.7336, LR test value: 23.353, p-value: 0.000001348
## Asymptotic standard error: 0.11018
## z-value: 6.6584, p-value: 0.000000000027681
## Wald statistic: 44.334, p-value: 0.000000000027681
##
## Log likelihood: -3439.88 for error model
## ML residual variance (sigma squared): 626.55, (sigma: 25.031)
## Number of observations: 741
## Number of parameters estimated: 5
## AIC: 6889.8, (AIC for lm: 6911.1)# spatial lag regression
# Note that the LM tests suggested lag model was a better choice here
fire_spatial_lag <- lagsarlm(FireDays ~ Calculated + Lag(FireDays,1), data=fire_model_data, W_dist_rowstd, zero.policy=T)
summary(fire_spatial_lag)##
## Call:lagsarlm(formula = FireDays ~ Calculated + Lag(FireDays, 1),
## data = fire_model_data, listw = W_dist_rowstd, zero.policy = T)
##
## Residuals:
## Min 1Q Median 3Q Max
## -95.2648 -13.8359 -9.5096 7.1668 130.9648
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.60646413 2.26212048 0.7102 0.477606
## Calculated 0.00493977 0.00039894 12.3822 < 0.00000000000000022
## Lag(FireDays, 1) 0.08661095 0.03240254 2.6730 0.007518
##
## Rho: 0.75736, LR test value: 30.342, p-value: 0.000000036214
## Asymptotic standard error: 0.096276
## z-value: 7.8665, p-value: 0.0000000000000035527
## Wald statistic: 61.882, p-value: 0.0000000000000036637
##
## Log likelihood: -3436.386 for lag model
## ML residual variance (sigma squared): 620.27, (sigma: 24.905)
## Number of observations: 741
## Number of parameters estimated: 5
## AIC: 6882.8, (AIC for lm: 6911.1)
## LM test for residual autocorrelation
## test value: 0.74624, p-value: 0.38767stargazer(olstime, fixed_effects_fire, fire_spatial_error, digits=1, type = "html", out="ols_ml_sp.html")##
## <table style="text-align:center"><tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="3"><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="3" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td colspan="3">FireDays</td></tr>
## <tr><td style="text-align:left"></td><td><em>OLS</em></td><td><em>panel</em></td><td><em>spatial</em></td></tr>
## <tr><td style="text-align:left"></td><td><em></em></td><td><em>linear</em></td><td><em>error</em></td></tr>
## <tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td><td>(3)</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Calculated</td><td>0.01<sup>***</sup></td><td>0.01<sup>***</sup></td><td>0.005<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.000)</td><td>(0.000)</td><td>(0.000)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Lag(FireDays, 1)</td><td>0.1<sup>***</sup></td><td></td><td>0.1<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.03)</td><td></td><td>(0.03)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">Constant</td><td>15.8<sup>***</sup></td><td></td><td>17.5<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(1.2)</td><td></td><td>(3.5)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>741</td><td>742</td><td>741</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.2</td><td>0.2</td><td></td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.2</td><td>0.2</td><td></td></tr>
## <tr><td style="text-align:left">Log Likelihood</td><td></td><td></td><td>-3,439.9</td></tr>
## <tr><td style="text-align:left">sigma<sup>2</sup></td><td></td><td></td><td>626.6</td></tr>
## <tr><td style="text-align:left">Akaike Inf. Crit.</td><td></td><td></td><td>6,889.8</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>25.6 (df = 738)</td><td></td><td></td></tr>
## <tr><td style="text-align:left">F Statistic</td><td>84.1<sup>***</sup> (df = 2; 738)</td><td>155.8<sup>***</sup> (df = 1; 737)</td><td></td></tr>
## <tr><td style="text-align:left">Wald Test</td><td></td><td></td><td>44.3<sup>***</sup> (df = 1)</td></tr>
## <tr><td style="text-align:left">LR Test</td><td></td><td></td><td>23.4<sup>***</sup> (df = 1)</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="3" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>