Ev2007Higod4

From InterSciWiki
Jump to: navigation, search

Printable Moral Gods Model & Data

Princomp variables were not tested. All but dx$SuperjhWriting are significant. Rsq=.40-.43 range in different runs.

Save Variable Names Wanted in gg

dx$v2007 FxCmtyWages AnimXbwealth No_rain_Dry SuperjhWriting v149
HiGod4      <-data.frame(dx$v2007,dx$FxCmtyWages,dx$AnimXbwealth,dx$No_rain_Dry,dx$SuperjhWriting,dx$v149)  #MAKE DATAFRAME
        HiGod4A    <-data.frame(dx$v2007,dx$FxCmtyWages,dx$AnimXbwealth,dx$No_rain_Dry,dx$SuperjhWriting,dx$v149)  #MAKE DATAFRAME
save(HiGod4A,file="HiGod4A.Rdata")   #http://facweb.knowlton.ohio-state.edu/pviton/courses2/crp87105/data-frame.html
dx$SuperjhWriting[172]=1         #modal 
dx$SuperjhWriting[115]=1         #modal 
        HiGod4Modal<-data.frame(dx$v2007,dx$FxCmtyWages,dx$AnimXbwealth,dx$No_rain_Dry,dx$SuperjhWriting,dx$v149)  #MAKE DATAFRAME
save(HiGod4Modal,file="HiGod4AModal.Rdata")   #http://facweb.knowlton.ohio-state.edu/pviton/courses2/crp87105/data-frame.html
dx$SuperjhWriting[172]=2         #median
dx$SuperjhWriting[115]=2         #median
        HiGod4Median<-data.frame(dx$v2007,dx$FxCmtyWages,dx$AnimXbwealth,dx$No_rain_Dry,dx$SuperjhWriting,dx$v149)  #MAKE DATAFRAME
save(HiGod4Median,file="HiGod4AMedian.Rdata")   #http://facweb.knowlton.ohio-state.edu/pviton/courses2/crp87105/data-frame.html
table(gg$v149,dx$v149)
table(gg$AnimXbwealth,dx$AnimXbwealth)
table(gg$No_rain_Dry,dx$No_rain_Dry)
table(gg$FxCmtyWages,dx$FxCmtyWages)
table(gg$SuperjhWriting,dx$SuperjhWriting)
table(gg$FxCmtyWages,dx$FxCmtyWages)  #n=184 the only table with NAs
table(gg$v149,dx$v149) #DIFFER   
    1  2  3  4  5
 1 33 17  7  2 14
 2 18 14  6  6  5
 3  5  3  4  1  8
 4  6  5  1  0  0
 5 11 10  3  3  4
table(gg$AnimXbwealth,dx$AnimXbwealth) #DIFFER
    0  1  2  3  4  5  7  8  9
 0 85  9 11 11  2  1  1  4  4
 1 12  1  1  1  1  0  0  0  0
 2 12  2  1  0  0  0  0  0  0
 3 10  2  0  1  0  0  0  0  0
 4  2  1  1  0  0  0  0  0  0
 5  1  0  0  0  0  0  0  0  0
 7  0  0  0  0  1  0  0  0  0
 8  3  1  0  0  0  0  0  0  0
 9  3  0  1  0  0  0  0  0  0
table(gg$No_rain_Dry,dx$No_rain_Dry) #DIFFER
    1  2  3  4
 1 54 19 22  7
 2 22  4  5  4
 3 17  9  4  3
 4  9  3  2  2
table(gg$FxCmtyWages,dx$FxCmtyWages) #SAME  n=184
     0   1
 0 161  12
 1  12   1
table(gg$SuperjhWriting,dx$SuperjhWriting) #DIFFER
     1  2  3  4  5  6  8 10
 1  31 18  5  4  1  5  5  3
 2  13 11  4  5  0  3  5  4
 3   4  2  1  3  0  0  1  0
 4   8  3  0  3  0  1  1  0
 5   0  0  0  0  0  0  0  1
 6   6  2  0  1  0  0  2  1
 8   7  2  2  1  0  2  1  1
table(gg$v149,dx$v149)  #DIFFER
    1  2  3  4  5
 1 33 17  7  2 14
 2 18 14  6  6  5
 3  5  3  4  1  8
 4  6  5  1  0  0
 5 11 10  3  3  4

using h$data[,"x"]

v2007MoralGods

Ev2007Higod - HiGodNoWy Bayesian network parameters for testing White et al vs. Eff and Brown Comments149 BE CAREFUL ABOUT THE SUBHEADING BELOW Ev2007Higod4#for_Bayesian_Graph Does this have all the Brown and Eff variables that work?

h[7]  Douglas R. White (talk) 06:53, 12 February 2014 (PST)
$OtherStats
  d l e Weak.Identification.Fstat R2.final.model R2.UR.model nimp nobs BClambda
1 1 0 0                  21.34552      0.4209838   0.4657769    2  186     none
h[4]
$Rmodel
                 coef  stdcoef     VIF  relimp    pval  hcpval star                                            desc
(Intercept)   0.14702      NaN     NaN     NaN 0.56803 0.54191                                                 <NA>
AnimXbwealth  0.10524  0.17851 1.30363 0.07216 0.00599 0.01821  ***                                    AnimXbwealth
FxCmtyWages   0.81479  0.17729 1.01409 0.03262 0.00197 0.00000  ***                                     FxCmtyWages
No_rain_Dry   0.16355  0.14167 1.44789 0.06918 0.03844 0.08866   **                                     No_rain_Dry
v149          0.16029  0.20017 1.38606 0.08445 0.00288 0.00621  ***                    Scale 1- Writing and Records
v1650        -0.01756 -0.10002 1.09125 0.00542 0.09829 0.09407    * Frequency of External Warfare (Resolved Rating)
Wy            0.63593  0.33636 1.54355 0.15992 0.00000 0.00001  ***                                Network lag term
library(mice)
library(foreign)
library(stringr)
library(AER)
library(spdep)
library(psych)
library(geosphere)
library(relaimpo)
library(linprog)
library(dismo)
library(forward)
library(pastecs)
library(classInt)
library(maps)
library(dismo)
library(plyr)
library(aod)
library(reshape)
 library(RColorBrewer)
 library(XML)
 library(tm)
 library(mlogit)
#The Dow-Eff functions, as well as the four ethnological datasets, are contained in an R-workspace, located in the cloud.
#load(url("http://dl.dropbox.com/u/9256203/DEf01.Rdata"), .GlobalEnv)
 load(url("http://dl.dropbox.com/u/9256203/DEf01d.Rdata"), .GlobalEnv)
 #load(url("http://dl.dropbox.com/u/9256203/DEf01f.Rdata"), .GlobalEnv)
 load(url("http://capone.mtsu.edu/eaeff/downloads/DEf01f.Rdata")) # uploaded Douglas R. White (talk) 08:20, 3 August 2014 (PDT)
ls()  #-can see the objects contained in DEf01.Rdata
#The setDS( xx ) command sets one of the four ethnological datasets as the source for the subsequent analysis. The four valid options for xx are: “WNAI”, “LRB”, “EA”, “SCCS”. The setDS() command creates objects:
setDS("SCCS")
dx$HiGod=dx$v2007
addesc("HiGod", "HiGod")  # originsymbol
dx$v2006v149=dx$2006*dx$v149
dx$FxCmtyWages=((dx$v2125==3)*1)*(dx$v61==6)*1 #v2125 is (new) Wages - v61 is Fixity of community
addesc("FxCmtyWages","FxCmtyWages") 
dx$SuperjhWriting<-dx$v237*(1+((dx$v149>=3)*1))
addesc("SuperjhWriting","SuperjhWriting") 
dx$AnimXbwealth=dx$v206*(dx$v208==1)*1
addesc("AnimXbwealth","AnimXbwealth") 
dx$No_rain_Dry<-dx$v2003+(dx$v855>4)*1
addesc("No_rain_Dry","No_rain_Dry")
#dx$sqHiGod<-dx$HiGod^2
#addesc("sqAXB","sqAXB")
dx$sq149<-dx$v149^2
addesc("sq149","sq149")
dx$sqAXB<-dx$AnimXbwealt^2
addesc("sqAXB","sqAXB")
dx$sqNRD<-dx$No_rain_Dry^2
addesc("sqNRD","sqNRD")
#mkdummy("v244", 7) ## [1] "This dummy variable is named v244.d7"
#mkdummy("v245", 2) ## [1] "This dummy variable is named v245.d2"  milk
#mkdummy("v233", 6) ## [1] "This dummy variable is named v233.d6"


#evm<-unique(c(polyg,"v279.d1","rectang",femecon,path,avoid)) List_of_Rgui_to_CoSSci_Models#Eff.27s_Oct_18_2013_MiVersion_of_DEf01cmkscale  Ev676.4
##evm<-unique(c(HiGod, "v2006v149","bio.5","HiGod","FxCmtyWages","SuperjhWriting","AnimXbwealth","No_rain_Dry","v149","v2003","sq149","sqAXB","sqNRD", "v2004","v2005","v2006","v272","v1683",                "v206",  "v1650","v1695","v270","v203","v204","v205","v21","v53","v665","v670","v7", "v872", "v157",   "PCAP","PCsize"))
##evm<-unique(c(HiGod, "bio.5",        "FxCmtyWages","SuperjhWriting","AnimXbwealth","No_rain_Dry","v149","v2003","sq149","sqAXB","sqNRD", "v2004","v2005","v2006","v272","v1683",                "v206",  "v1650","v1695","v270","v203","v204","v205","v21","v53","v665","v670","v7", "v872", "v157",   "PCAP","PCsize"))
##Error in unique(c(HiGod, "bio.5", "HiGod", "FxCmtyWages", "SuperjhWriting",  :  error in evaluating the argument 'x' in selecting a method for function 'unique': Error: object 'PCAP' not found
#evm <- c("PCAP","PCsize","bio.5","HiGod","FxCmtyWages","SuperjhWriting","AnimXbwealth","No_rain_Dry","v149","v2003","sq149","sqAXB","sqNRD", "v2004","v2005","v2006","v272","v1683",                "v206",  "v1650","v1695","v270","v203","v204","v205","v21","v53","v665","v670","v7", "v872", "v157")            
#ABOVE AND BELOW RUN BUT DONT PRODUCE PCAP or PCsize
evm <- c("HiGod","PCAP","PCsize","bio.5","FxCmtyWages","SuperjhWriting","AnimXbwealth","No_rain_Dry","v149","v2003","sq149","sqAXB","sqNRD", "v2004","v2005","v2006","v272","v1683",                "v206",  "v1650","v1695","v270","v203","v204","v205","v21","v53","v665","v670","v7", "v872", "v157")  
# Brown Eff # fec$scales  
# type="pc1", add.descrip="1st PC: Agricultural potential high")
# type="pc1", add.descrip="1st PC: Community size large")
#v1695 foodscarc Chronic resource scarcity high  BROWN EFF FULLY INCLUDED EXCEPT FOR 1st PCAs
#v206 anim Percentage subsistence: Animal husbandry
#v270 caststrat Degree of case stratification
#v1650 eeextwar Frequency of external war high
smi <- doMI(evm, nimp = 2, maxit = 3)
#==agricultural potential delete if the above script runs==
# set type="pc1" for first principal component. To ensure that the variable has the 
# expected meaning, use the "set.direction" option.
 fec<-mkscale(compvarbs="agp", udnavn="PCAP", impdata=smi, set.direction="v921",
             type="pc1", add.descrip="1st PC: Agricultural potential high")
#--check reasonableness of scale--
fec$stats
fec$corrs
smi[,names(fec$scales)]<-fec$scales
# ==commmunity size==
fec<-mkscale(compvarbs="csz", udnavn="PCsize", impdata=smi, set.direction="v63",
             type="pc1", add.descrip="1st PC: Community size large")
#--check reasonableness of scale--
fec$stats
fec$corrs
smi[,names(fec$scales)]<-fec$scales
### END ==ag
# --dependent variable--
dpV <- "HiGod" 
# --independent variables in UNrestricted model--
UiV <- c("PCAP","PCsize","v1650","v1695","v270","bio.5","FxCmtyWages","SuperjhWriting","AnimXbwealth","No_rain_Dry","sqAXB","sqNRD","v149","v2003","v272","v1683","v2004", "v2005",       "v203", "v204","v205","v21","v53","v665","v670", "v7")  #,"v2125"
# --additional exogenous variables (use in Hausman tests)--
oxog <- c("v2006","v149","v2006v149")
# --independent variables in restricted model (all must be in UiV above)--
RiV <- c("v2006v149","PCAP","PCsize","v1650","v270","FxCmtyWages","AnimXbwealth","No_rain_Dry")#,"v149")     #,"v272","v1650","v1683","sqAXB","sqNRD","SuperjhWriting" NS ,"v2125","v2004","v2005","v670","v1695" 
#WHY AnimXbwealth	FxCmtyWages	No_rain_Dry	v149	v1650	v270
h <- doOLS(smi, depvar = dpV, indpv = UiV, rindpv = RiV, othexog = oxog, dw = TRUE,  lw = TRUE, ew = FALSE, stepW = TRUE, relimp = TRUE, slmtests = FALSE, full.set=FALSE)
#Anthon: Douglas R. White (talk) 12:50, 3 February 2014 (PST) "full.set=FALSE"
CSVwrite(h[1-11], "Ev2007Higod4.2125", FALSE)
#v1695 foodscarc Chronic resource scarcity high  BROWN EFF FULLY INCLUDED EXCEPT FOR 1st PCAs
#v206 anim Percentage subsistence: Animal husbandry  SKIP
#v270 caststrat Degree of case stratification
#v1650 eeextwar Frequency of external war high

for Bayesian Graph

V2181_Reincarnation#for_Bayesian_Graph
#WNAI4AvoidWiMo
Wy <- h$data[,"Wy"]  #this gives Wy
HiGod <- h$data[,"HiGod"] 
FxCmtyWages <- h$data[,"FxCmtyWages"] 
AnimXbwealth <- h$data[,"AnimXbwealth"] 
No_rain_Dry <- h$data[,"No_rain_Dry"] 
Writing <- h$data[,"v149"] 
PCAP <- h$data[,"PCAP"] 
PCsize <- h$data[,"PCsize"] 
v1695 <- h$data[,"v1695"] 
v270 <- h$data[,"v270"] 
v1650 <- h$data[,"v1650"] 
out <- lm(HiGod ~ Wy + FxCmtyWages + AnimXbwealth + No_rain_Dry + v149 )
summary(out)
out2 <- lm(HiGod ~ FxCmtyWages + AnimXbwealth + No_rain_Dry + v149 )
summary(out)
setwd("/Users/drwhite/")
HHiGGod <- data.frame(Wy,HiGod,FxCmtyWages,AnimXbwealth,No_rain_Dry,Writing,v1695,v270,v1650)             #,PCAP,PCsizesave(HHiGGod,file="HHiGGod.Rda")
#That was at setwd("/Users/drwhite/desktop")
#Other stuff was at setwd("/Users/drwhite/")
write.csv(HHiGGod,file="HHiGGod.csv")
 load('HiGod4csvR.Rdata')
   h2=HHiGGod         #HiGod4csvR
#      X        Wy HiGod FxCmtyWages AnimXbwealth No_rain_Dry v1695 v270 v1650
#1     1 2.3319351     1           0            0           4   2.0    2  17.0
#2     2 2.1613937     3           0            0           3   3.0    1   1.0
#getwd()
#[1] "/Users/drwhite"
library(bnlearn)
# setwd("/Users/drwhite/desktop")
load('HiGod4csvR2.Rdata')
   h2=HHiGGod    
h2[]=lapply(h2,factor)         #tried run with no factors just real numbers but that didn't work either
str(h2)
res = gs(h2)
plot(res)
#save.csv("Ev2007.csv",file=
#save.csv("LRBsortPlusWvRecod3.csv",file="

Bayesian Network Working Version

Ev2007Higod4#Bayesian Network Working Version

setwd("/Users/drwhite/")
HiGodNoWy <- data.frame(HiGod,FxCmtyWages,AnimXbwealth,No_rain_Dry,Writing,v1695,v270,v1650)             #,PCAP,PCsizesave(HHiGGod,file="HHiGGod.Rda")
#That was at setwd("/Users/drwhite/desktop")
#Other stuff was at setwd("/Users/drwhite/")
write.csv(HiGodNoWy,file="HiGodNoWy.csv")
 load('HiGodNoWycsvR.Rdata')
   h2=HiGodNoWy         #HiGod4csvR
#      X        Wy HiGod FxCmtyWages AnimXbwealth No_rain_Dry v1695 v270 v1650
#1     1 2.3319351     1           0            0           4   2.0    2  17.0
#2     2 2.1613937     3           0            0           3   3.0    1   1.0
#getwd()
#[1] "/Users/drwhite"
library(bnlearn)
# setwd("/Users/drwhite/desktop")
load('HiGodNoWy.Rdata')
   h2=HiGodNoWy    
h2[]=lapply(h2,factor)         #tried run with no factors just real numbers but that didn't work either
str(h2)
res = gs(h2)
plot(res)
bn.gs <- gs(HiGodNoWy)            ##NO ERROR after eliminating Wy
res2   = pdag2dag(bn.gs, ordering=c("FxCmtyWages","AnimXbwealth","HiGod","Writing","No_rain_Dry",'v1695','v270','v1650'))
res2
fit1=bn.fit(res2,HiGodNoWy)
fit1

SEE: HiGodNoWy Bayesian network parameters

CSVwrite(h, "Ev2007Higod4_d", FALSE)  #totry  "bio.2","v2006","v53","bio.5","v1650","v272","v53","sqAXB","HiGodSq","

Result

$DependVarb "Dependent variable='HiGod': HiGod"
RiV <- c("FxCmtyWages","AnimXbwealth","No_rain_Dry","v149")  
$didwell "FxCmtyWages"  "No_rain_Dry"  "v149"         "AnimXbwealth"
> h[5]
$Diagnostics
                                                         Fstat           df pvalue star
RESET test. H0: model has correct functional form        3.379 1.222560e+16  0.066    * <- modest problem: sq corrections failed
Wald test. H0: appropriate variables dropped             0.000 6.165094e+06  1.000     
Breusch-Pagan test. H0: residuals homoskedastic          7.696 2.349440e+07  0.006  ***
Shapiro-Wilkes test. H0: residuals normal                4.528 1.652762e+07  0.033   **
Hausman test. H0: Wy is exogenous                        0.000 1.862600e+04  1.000     
Sargan test. H0: residuals uncorrelated with instruments 0.896 5.268058e+07  0.344     
later> h[3] $Rmodel
                 coef stdcoef     VIF  relimp  hcpval    pval star desc
(Intercept)  -0.36382     NaN     NaN     NaN 0.20729 0.23165      <NA>
AnimXbwealth  0.09997 0.16958 1.30990 0.07071 0.02500 0.01031   ** AnimXbwealth
FxCmtyWages   0.81484 0.17730 1.01236 0.03244 0.00000 0.00228  *** FxCmtyWages
No_rain_Dry   0.15425 0.13362 1.42521 0.06858 0.11944 0.05260    * No_rain_Dry
v149          0.14687 0.18341 1.22750 0.08356 0.00952 0.00415  *** Writing and Records-Scale 1
Wy            0.82526 0.36401 1.43621 0.17774 0.00000 0.00000  *** Network lag term
earlier > h[3]
$Rmodel
                 coef stdcoef     VIF  relimp  hcpval    pval star                         desc
(Intercept)  -0.42434     NaN     NaN     NaN 0.13609 0.17249                              <NA> EARLIER MODEL
AnimXbwealth  0.11068 0.18775 1.26441 0.07602 0.01089 0.00395  *** AnimXbwealth
FxCmtyWages   0.83570 0.18184 1.01247 0.03333 0.00000 0.00181  *** FxCmtyWages
v149          0.15506 0.19364 1.21220 0.08803 0.00621 0.00240  *** Scale 1- Writing and Records
v2003         0.19695 0.11057 1.30849 0.05287 0.16813 0.09520    *   Rain
Wy            0.83115 0.36874 1.37997 0.17921 0.00000 0.00000  *** Network lag term
> h[6]
$OtherStats
  d l e Weak.Identification.Fstat R2.final.model R2.UR.model nimp nobs
1 1 0 0                  15.07052      0.3959725   0.4187296    2  186
Optimal weight matrix was W=1.0*Wd+.00*Wl
> h[4]
$EndogeneityTests  Hausman
             weakidF p.Sargan n.IV Fstat         df pvalue star
AnimXbwealth 149.843    0.428    8 0.083 7.645830e+05  0.773     
FxCmtyWages    3.650    1.000    4 0.034 7.689262e+07  1.000     
No_rain_Dry   25.000    0.651   10 1.863 1.327025e+13  0.172     
v149           6.338    0.506   12 2.000 2.240400e+07  0.137

Codes from Brown and Eff

Brown_and_Eff_2010_code_-_Mac -Christian Brown *Christian Brown and Anthon Eff, 2010 pdf or The State and the Supernatural: Support for Prosocial Behavior, Structure and Dynamics: eJournal of Anthropological Sciences, 4(1) art 1. http://escholarship.org/uc/item/5rh6z6z6pdf NOT THE RIGHT URL! http://intersci.ss.uci.edu/wiki/pdf/Brown&EffTheStateandtheSupernatural-SupportforProsocialBehavior.pdf searchable text when downloaded

v921, v928 PCAP 1st PC: Agricultural potential high
v63, v237 PCsize 1st PC: Community size large   /// PCsizeSq PCsize squared   
v1695 foodscarc Chronic resource scarcity high
v206 anim Percentage subsistence: Animal husbandry
v270 caststrat Degree of case stratification
v1650 eeextwar Frequency of external war high

  • large set
AP1=SCCS$v921,AP2=SCCS$v928,
foodscarc=SCCS$v1685,
cultints=SCCS$v232,
anim=SCCS$v206,
milk=(SCCS$v245>1)*1,
classtrat=SCCS$v270, 
caststrat=SCCS$v272, 
superjh=SCCS$v237, 
money=SCCS$v155,
exogamy=SCCS$v72,
commland=SCCS$v1726,
inhreal=(SCCS$v278>1)*1,
inhmove=(SCCS$v279>1)*1,
marrgood=(SCCS$v208<4)*1,
commsize=SCCS$v63,
popdens=SCCS$v64,
homicide=SCCS$v1665, 
assault=SCCS$v1666, 
theft=SCCS$v1667,
persviol=SCCS$v666,
winconfl=SCCS$v767,
bwnconfl=SCCS$v768,
physforce=SCCS$v770,
intwar=SCCS$v773,
frqintwar=SCCS$v891,
eeintwar=SCCS$v1649, 
eeextwar=SCCS$v1650 )