Example 1
## Warning: package 'spsur' was built under R version 3.6.2
data("spc", package = "spsur")
head(spc)
## COUNTY WAGE83 UN83 NMR83 SMSA WAGE82 WAGE81 UN80
## 1 UNION 1.003127 0.080500 -0.002217 1 1.108662 1.146178 0.130375
## 2 DELAWARE 1.039972 0.122174 0.018268 1 1.071271 1.104241 0.189603
## 3 LICKING 1.050196 0.095821 -0.013681 1 1.058375 1.094732 0.124125
## 4 MADISON 1.052210 0.102941 -0.000959 1 1.049791 1.065032 0.139394
## 5 FRANKLIN 1.055406 0.111965 -0.005160 1 1.075883 1.085075 0.183311
## 6 FAIRFIELD 1.048299 0.095192 -0.008934 1 1.056021 1.110068 0.155076
## NMR80 WAGE80
## 1 -0.010875 1.084886
## 2 0.041886 1.110426
## 3 -0.004158 1.069776
## 4 -0.004271 1.067895
## 5 -0.015568 1.079768
## 6 0.004514 1.099969
spcformula <- WAGE83 | WAGE81 ~ UN83 + NMR83 + SMSA | UN80 + NMR80 + SMSA
LMs.spc <- lmtestspsur(Form = spcformula, data = spc, W = Wspc)
## LM-Stat. DF p-value
## LM-SUR-SLM 5.2472 2 0.0725 .
## LM-SUR-SEM 3.3050 2 0.1916
## LM*-SUR-SLM 2.1050 2 0.3491
## LM*-SUR-SEM 0.1628 2 0.9218
## LM-SUR-SARAR 5.7703 4 0.2170
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
spcsur.slm <- spsurml(Form = spcformula, data = spc, type = 'slm', W = Wspc)
## Initial point: log_lik: 113.197 lambdas: -0.472 -0.446
## Iteration: 1 log_lik: 114.085 lambdas: -0.506 -0.482
## Iteration: 2 log_lik: 114.096 lambdas: -0.506 -0.482
## Time to fit the model: 2.25 seconds
## Computing marginal test...
## Time to compute covariances: 0.39 seconds
## Call:
## spsurml(Form = spcformula, data = spc, W = Wspc, type = "slm")
##
##
## Spatial SUR model type: slm
##
## Equation 1
## Estimate Std. Error t value Pr(>|t|)
## (Intercept)_1 1.4955217 0.2467240 6.0615 5.183e-07 ***
## UN83_1 0.8070029 0.2557439 3.1555 0.003179 **
## NMR83_1 -0.5194114 0.2590550 -2.0050 0.052318 .
## SMSA_1 -0.0073247 0.0118519 -0.6180 0.540347
## lambda_1 -0.5057334 0.2405734 -2.1022 0.042401 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## R-squared: 0.6224
## Equation 2
## Estimate Std. Error t value Pr(>|t|)
## (Intercept)_2 1.7094414 0.2925620 5.8430 1.024e-06 ***
## UN80_2 -0.6745562 0.3870737 -1.7427 0.08969 .
## NMR80_2 0.7502934 0.3842670 1.9525 0.05847 .
## SMSA_2 0.0014181 0.0241859 0.0586 0.95356
## lambda_2 -0.4821428 0.2557758 -1.8850 0.06730 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## R-squared: 0.4743
## Variance-Covariance Matrix of inter-equation residuals:
## 0.0003085954 -0.0003561928
## -0.0003561928 0.0015864976
## Correlation Matrix of inter-equation residuals:
## 1.000000 -0.509062
## -0.509062 1.000000
##
## R-sq. pooled: 0.6603
## Log-Likelihood: 114.096
## Breusch-Pagan: 6.516 p-value: (0.0107)
## LMM: 0.50489 p-value: (0.477)
spcSUR.sim <- spsurml(Form = spcformula, data = spc, type = 'sim', W = Wspc)
## Initial point:
## log_lik: 110.423
## Iteration: 1 log_lik: 111.201
## Iteration: 2 log_lik: 111.348
## Iteration: 3 log_lik: 111.378
## Time to fit the model: 0.11 seconds
## Time to compute covariances: 0.02 seconds
spcSUR.slx <- spsurml(Form = spcformula, data = spc, type = 'slx', W = Wspc)
## Initial point:
## log_lik: 113.659
## Iteration: 1 log_lik: 115.304
## Iteration: 2 log_lik: 115.792
## Iteration: 3 log_lik: 115.915
## Iteration: 4 log_lik: 115.942
## Time to fit the model: 0.15 seconds
## Time to compute covariances: 0.01 seconds
spcSUR.sem <- spsurml(Form = spcformula, data = spc, type = 'sem', W = Wspc)
## Initial point: log_lik: 112.821 deltas: -0.556 -0.477
## Iteration: 1 log_lik: 113.695 rhos: -0.618 -0.537
## Iteration: 2 log_lik: 113.719 rhos: -0.628 -0.548
## Time to fit the model: 3.03 seconds
## Computing marginal test...
## Time to compute covariances: 0.36 seconds
spcSUR.sarar <- spsurml(Form = spcformula, data = spc, type = 'sarar', W = Wspc)
## Initial point: log_lik: 113.406 lambdas: -0.343 -0.526 rhos: -0.28 0.113
## Iteration: 1 log_lik: 114.574 lambdas: -0.384 -0.783 rhos: -0.307 0.411
## Iteration: 2 log_lik: 114.824 lambdas: -0.394 -0.905 rhos: -0.303 0.553
## Iteration: 3 log_lik: 115.019 lambdas: -0.395 -0.995 rhos: -0.296 0.657
## Iteration: 4 log_lik: 115.187 lambdas: -0.388 -1 rhos: -0.271 0.686
## Iteration: 5 log_lik: 115.237 lambdas: -0.381 -1 rhos: -0.254 0.699
## Iteration: 6 log_lik: 115.26 lambdas: -0.376 -1 rhos: -0.243 0.707
## Time to fit the model: 34.7 seconds
## Time to compute covariances: 0.2 seconds
spcSUR.sdm <- spsurml(Form = spcformula, data = spc, type = 'sdm', W = Wspc)
## Initial point: log_lik: 115.602 lambdas: -0.496 -0.454
## Iteration: 1 log_lik: 116.566 lambdas: -0.512 -0.472
## Iteration: 2 log_lik: 116.646 lambdas: -0.482 -0.449
## Iteration: 3 log_lik: 116.695 lambdas: -0.454 -0.428
## Time to fit the model: 2.33 seconds
## Computing marginal test...
## Time to compute covariances: 0.15 seconds
spcSUR.sdem <- spsurml(Form = spcformula, data = spc, type = 'sdem', W = Wspc)
## Initial point: log_lik: 116.323 deltas: -0.625 -0.695
## Iteration: 1 log_lik: 118.191 rhos: -0.745 -0.925
## Iteration: 2 log_lik: 118.442 rhos: -0.761 -1
## Iteration: 3 log_lik: 118.489 rhos: -0.752 -1
## Time to fit the model: 3.11 seconds
## Computing marginal test...
## Time to compute covariances: 0.2 seconds
spcSUR.slm.3sls <-spsur3sls(Form = spcformula, data = spc, type = 'slm', W = Wspc)
## Time to fit the model: 0.03 seconds
#> Time to fit the model: 0.17 seconds
summary(spcSUR.slm.3sls)
## Call:
## spsur3sls(Form = spcformula, data = spc, W = Wspc, type = "slm")
##
##
## Spatial SUR model type: slm
##
## Equation 1
## Estimate Std. Error t value Pr(>|t|)
## (Intercept)_1 1.4105098 0.3395643 4.1539 0.0001849 ***
## UN83_1 0.8271804 0.2777129 2.9785 0.0050887 **
## NMR83_1 -0.5604625 0.2915888 -1.9221 0.0623193 .
## SMSA_1 -0.0096794 0.0133180 -0.7268 0.4719261
## lambda_1 -0.4242869 0.3314872 -1.2799 0.2085301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## R-squared: 0.6119
## Equation 2
## Estimate Std. Error t value Pr(>|t|)
## (Intercept)_2 1.7217295 0.3866951 4.4524 7.544e-05 ***
## UN80_2 -0.6807387 0.4270516 -1.5940 0.11943
## NMR80_2 0.8270205 0.4294552 1.9257 0.06185 .
## SMSA_2 0.0026117 0.0276672 0.0944 0.92530
## lambda_2 -0.4936829 0.3397513 -1.4531 0.15463
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## R-squared: 0.4767
## Variance-Covariance Matrix of inter-equation residuals:
## 0.0003567306 -0.0003936538
## -0.0003936538 0.0019083726
## Correlation Matrix of inter-equation residuals:
## 1.0000000 -0.4771036
## -0.4771036 1.0000000
##
## R-sq. pooled: 0.6628
Example 2
library(foreign)
library(systemfit)
## Warning: package 'systemfit' was built under R version 3.6.2
## Loading required package: Matrix
## Loading required package: car
## Warning: package 'car' was built under R version 3.6.2
## Loading required package: carData
## Warning: package 'carData' was built under R version 3.6.1
## Loading required package: lmtest
## Warning: package 'lmtest' was built under R version 3.6.2
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 3.6.2
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Please cite the 'systemfit' package as:
## Arne Henningsen and Jeff D. Hamann (2007). systemfit: A Package for Estimating Systems of Simultaneous Equations in R. Journal of Statistical Software 23(4), 1-40. http://www.jstatsoft.org/v23/i04/.
##
## If you have questions, suggestions, or comments regarding the 'systemfit' package, please use a forum or 'tracker' at systemfit's R-Forge site:
## https://r-forge.r-project.org/projects/systemfit/
hsb2 <- read.dta("https://stats.idre.ucla.edu/stat/stata/notes/hsb2.dta")
head(hsb2)
## id female race ses schtyp prog read write math science socst
## 1 70 male white low public general 57 52 41 47 57
## 2 121 female white middle public vocation 68 59 53 63 61
## 3 86 male white high public general 44 33 54 58 31
## 4 141 male white high public vocation 63 44 47 53 56
## 5 172 male white middle public academic 47 52 57 53 61
## 6 113 male white middle public academic 44 52 51 63 61
## id female race ses schtyp
## Min. : 1.00 male : 91 hispanic : 24 low :47 public :168
## 1st Qu.: 50.75 female:109 asian : 11 middle:95 private: 32
## Median :100.50 african-amer: 20 high :58
## Mean :100.50 white :145
## 3rd Qu.:150.25
## Max. :200.00
## prog read write math science
## general : 45 Min. :28.00 Min. :31.00 Min. :33.00 Min. :26.00
## academic:105 1st Qu.:44.00 1st Qu.:45.75 1st Qu.:45.00 1st Qu.:44.00
## vocation: 50 Median :50.00 Median :54.00 Median :52.00 Median :53.00
## Mean :52.23 Mean :52.77 Mean :52.65 Mean :51.85
## 3rd Qu.:60.00 3rd Qu.:60.00 3rd Qu.:59.00 3rd Qu.:58.00
## Max. :76.00 Max. :67.00 Max. :75.00 Max. :74.00
## socst
## Min. :26.00
## 1st Qu.:46.00
## Median :52.00
## Mean :52.41
## 3rd Qu.:61.00
## Max. :71.00
r1 <- read~female + as.numeric(ses) + socst
r2 <- math~female + as.numeric(ses) + science
fitsur <- systemfit(list(readreg = r1, mathreg = r2), data=hsb2)
summary(fitsur)
##
## systemfit results
## method: OLS
##
## N DF SSR detRCov OLS-R2 McElroy-R2
## system 400 392 22835.2 3227.86 0.405103 0.342707
##
## N DF SSR MSE RMSE R2 Adj R2
## readreg 200 196 12550.9 64.0351 8.0022 0.400037 0.390854
## mathreg 200 196 10284.4 52.4712 7.2437 0.411171 0.402158
##
## The covariance matrix of the residuals
## readreg mathreg
## readreg 64.0351 11.4952
## mathreg 11.4952 52.4712
##
## The correlations of the residuals
## readreg mathreg
## readreg 1.00000 0.19831
## mathreg 0.19831 1.00000
##
##
## OLS estimates for 'readreg' (equation 1)
## Model Formula: read ~ female + as.numeric(ses) + socst
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.6824980 2.9789550 6.94287 5.5019e-11 ***
## femalefemale -1.5111280 1.1510793 -1.31279 0.19079
## as.numeric(ses) 1.2183658 0.8399004 1.45061 0.14849
## socst 0.5699327 0.0562967 10.12373 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.002195 on 196 degrees of freedom
## Number of observations: 200 Degrees of Freedom: 196
## SSR: 12550.883066 MSE: 64.035118 Root MSE: 8.002195
## Multiple R-Squared: 0.400037 Adjusted R-Squared: 0.390854
##
##
## OLS estimates for 'mathreg' (equation 2)
## Model Formula: math ~ female + as.numeric(ses) + science
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.305181 2.998047 6.43925 9.0557e-10 ***
## femalefemale 1.160903 1.041641 1.11449 0.266432
## as.numeric(ses) 1.399639 0.742390 1.88531 0.060867 .
## science 0.575330 0.054328 10.58993 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.243704 on 196 degrees of freedom
## Number of observations: 200 Degrees of Freedom: 196
## SSR: 10284.364144 MSE: 52.471246 Root MSE: 7.243704
## Multiple R-Squared: 0.411171 Adjusted R-Squared: 0.402158
library(car)
restriction <- "readreg_femalefemale- mathreg_femalefemale"
linearHypothesis(fitsur, restriction, test = "Chisq")
## Linear hypothesis test (Chi^2 statistic of a Wald test)
##
## Hypothesis:
## readreg_femalefemale - mathreg_femalefemale = 0
##
## Model 1: restricted model
## Model 2: fitsur
##
## Res.Df Df Chisq Pr(>Chisq)
## 1 393
## 2 392 1 2.9626 0.08521 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Example 3
setwd("C:/Users/subas/Syncplicity/MyProjects_IMP/MY_Papers_V2/TRB 2021/EScotter_BayesianRule/")
it01 <- read.csv("IT_aadtMaster.csv")
mn= subset(it01, State=="MN")
dim(mn)
## [1] 11498 86
mn01 <- mn[, c("Default_AADT", "HU", "Pop" ,"C_Pop", "C_HU" , "Pop_Empl", "WAC", "RAC",
"Agg_Inc", "Agg_Veh", "Empl" )]
mn02= na.omit(mn01)
head(mn02)
## Default_AADT HU Pop C_Pop C_HU Pop_Empl WAC RAC Agg_Inc Agg_Veh Empl
## 30019 850 713 922 935 635 1544 770 486 48480900 690 622
## 30020 75 486 1270 1154 476 2044 1092 603 48077900 1154 774
## 30021 75 486 1270 1154 476 2044 1092 603 48077900 1154 774
## 30022 75 486 1270 1154 476 2044 1092 603 48077900 1154 774
## 30023 75 537 1269 1155 490 1979 240 666 47370100 1183 710
## 30024 75 537 1269 1155 490 1979 240 666 47370100 1183 710
r1 <- Pop~C_Pop+ C_HU+ Pop_Empl+Empl
r2 <- Default_AADT~Pop + HU + WAC+RAC+Agg_Veh+Agg_Inc
fitsur <- systemfit(list(readreg = r1, mathreg = r2), data=mn02)
summary(fitsur)
##
## systemfit results
## method: OLS
##
## N DF SSR detRCov OLS-R2 McElroy-R2
## system 21560 21548 536015438 0 0.950365 1
##
## N DF SSR MSE RMSE R2 Adj R2
## readreg 10780 10775 0 0.0 0.000 1.000000 1.000000
## mathreg 10780 10773 536015438 49755.4 223.059 0.496781 0.496501
##
## The covariance matrix of the residuals
## readreg mathreg
## readreg 3.22880e-24 -3.58957e-11
## mathreg -3.58957e-11 4.97554e+04
##
## The correlations of the residuals
## readreg mathreg
## readreg 1.0000000 -0.0902733
## mathreg -0.0902733 1.0000000
##
##
## OLS estimates for 'readreg' (equation 1)
## Model Formula: Pop ~ C_Pop + C_HU + Pop_Empl + Empl
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.95808e-13 3.89402e-14 -2.04367e+01 < 2.22e-16 ***
## C_Pop 4.88498e-15 9.81889e-17 4.97508e+01 < 2.22e-16 ***
## C_HU -5.89806e-16 9.12967e-17 -6.46032e+00 1.0899e-10 ***
## Pop_Empl 1.00000e+00 1.32414e-16 7.55204e+15 < 2.22e-16 ***
## Empl -1.00000e+00 2.94578e-16 -3.39469e+15 < 2.22e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0 on 10775 degrees of freedom
## Number of observations: 10780 Degrees of Freedom: 10775
## SSR: 0 MSE: 0 Root MSE: 0
## Multiple R-Squared: 1 Adjusted R-Squared: 1
##
##
## OLS estimates for 'mathreg' (equation 2)
## Model Formula: Default_AADT ~ Pop + HU + WAC + RAC + Agg_Veh + Agg_Inc
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.85640e+02 4.99827e+00 57.14772 < 2e-16 ***
## Pop 1.90870e-01 1.22159e-02 15.62470 < 2e-16 ***
## HU -9.24460e-03 1.13196e-02 -0.81669 0.41413
## WAC 2.96466e-02 1.27771e-03 23.20293 < 2e-16 ***
## RAC 6.36284e-01 1.80039e-02 35.34148 < 2e-16 ***
## Agg_Veh -7.98135e-01 1.23763e-02 -64.48897 < 2e-16 ***
## Agg_Inc 3.79471e-06 1.32418e-07 28.65692 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 223.059292 on 10773 degrees of freedom
## Number of observations: 10780 Degrees of Freedom: 10773
## SSR: 536015437.50521 MSE: 49755.447647 Root MSE: 223.059292
## Multiple R-Squared: 0.496781 Adjusted R-Squared: 0.496501