CreateModelDRWpolygyny.R

From InterSciWiki
Jump to: navigation, search
fixed lack of marrcaptives, plunder in indep_vars Doug 13:11, 19 January 2011 (PST)
fixed Whyte595=sccs$v594 --> Whyte595=sccs$v595,
Mac: setwd
setwd("/Users/Tara/Desktop/drwhite5files/sccs") #tara
setwd("/Users/drwhite/Documents/sccs") #this would work for the mac ignoring the next setwd, for your subdirectory on your mac
1A and 1B prototypes UPDATED Doug 12:16, 4 October 2010 (PDT)

1A CreateModelDRWpolygyny2.R

setwd("C:/My Documents/sccs")                 #PC
#setwd("/Users/drwhite/Documents/sccs") #Macbook <===================== comment out for PC!!!

library(sccs)
data(sccs)                                  #Run CreateModelDRWpolygyny2.R
depvar=sccs$v860
my_sccs<-data.frame(
dep_var=sccs$v860,
socname=sccs$socname,socID=sccs$"sccs#",
famsize=sccs$v68,    # similar to v80
exogamy=sccs$v72,
#famsizeB=sccs$v80,   # similar to v68
money=sccs$v155,
popdens=sccs$v156,
premarsexatt=sccs$v165,
premarsexfrq=sccs$v166,
malesexag=sccs$v175,
ndrymonth=sccs$v196,
gath=sccs$v203,
hunt=sccs$v204,
fish=sccs$v205,
anim=sccs$v206,
brideprice=(sccs$v208==1)*1,
nuclearfam=(sccs$v210<=3)*1,
ncmallow=sccs$v227,
cultints=sccs$v232,
tree=(sccs$v233==4)*1,
roots=(sccs$v233==5)*1,
cereals=(sccs$v233==6)*1,
settype=sccs$v234,
localjh=sccs$v236-1,
superjh=sccs$v237,
moralgods=sccs$v238,
segadlboys=sccs$v242,
plow=(sccs$v243>1)*1,
pigs=(sccs$v244==2)*1,
bovines=(sccs$v244==7)*1,
milk=(sccs$v245>1)*1,
agrlateboy=sccs$v300,
valchild=(sccs$v473+sccs$v474+sccs$v475+sccs$v476),
fratgrpstr=sccs$v570,
Whyte577=sccs$v577,    #take care using Whyte variables - only coded 1/2 the sample
Whyte580=sccs$v580,
Whyte584=sccs$v583,
Whyte585=sccs$v584,
Whyte595=sccs$v594,
Whyte602=sccs$v602,
Whyte615=sccs$v615,
Whyte620=sccs$v620,
Whyte626=sccs$v626,
Whyte629=sccs$v629,
Whyte630=sccs$v630,
Whyte631=sccs$v631,
Whyte632=sccs$v632,
Whyte633=sccs$v633,
Whyte635=sccs$v635,
#                    #take care using Paige variables - coded less than 1/2 the sample
Paige657=sccs$v657,  # summed in v663#  Paige657 Paige658 Paige659 Paige660 Paige661 Paige662
femproduceND=sccs$v658, #Paige658=sccs$v658,  # summed in v663 
Paige659=sccs$v659,  # summed in v663
Paige660=sccs$v660,  # summed in v663
Paige661=sccs$v661,  # summed in v663
Paige662=sccs$v662,  # summed in v663
fempower=sccs$v663, # sum of v657-662
interperviol=sccs$v666,   ###synonyms: violence / interviol / freintovio
migr=(sccs$v677==2)*1,
#                    #take care using Sanday variables - only coded 1/2 the sample
Sanday664=sccs$v664,  # summed in v669
Sanday665=sccs$v665,  # summed in v669
Sanday666=sccs$v666,  # summed in v669
Sanday667=sccs$v667,  # summed in v669
Sanday668=sccs$v668,  # summed in v669
Sanday669=sccs$v669, # sum of v664-668 
 #WHYTE Data Quality Whyte718 Whyte719 Whyte720 Whyte721 Whyte722 Whyte723 Whyte724 Whyte725
Whyte718=sccs$v718,    #take care using Whyte variables - only coded 1/2 the sample
Whyte719=sccs$v719,
Whyte720=sccs$v720,
Whyte721=sccs$v721,
Whyte722=sccs$v722,
Whyte723=sccs$v723,
Whyte724=sccs$v724,
Whyte725=sccs$v725,
 #Rohner Data Quality Codes 
Rohner798=sccs$v798,
Rohner799=sccs$v799,
Rohner800=sccs$v800,
Rohner801=sccs$v801,
Rohner802=sccs$v802,
Rohner803=sccs$v803,
Rohner804=sccs$v804,
Rohner805=sccs$v805,
Rohner806=sccs$v806,
Rohner807=sccs$v807,
Rohner808=sccs$v808,
Rohner809=sccs$v809,
Rohner810=sccs$v810,
Rohner811=sccs$v811,
Rohner812=sccs$v812,
Rohner813=sccs$v813,
foodtrade=sccs$v819,
fem_agri=sccs$v821, 
dateobs=sccs$v838,
rain=sccs$v854,
temp=sccs$v855,
#ecorich=sccs$v857,
ecorich=(sccs$v857==3|sccs$v857==4)*1+(sccs$v857==5)*2,
pctFemPolyg=sccs$v872,
marrcaptives=sccs$v870,
femsubs=sccs$v890,
intwar=sccs$v891,    # similar to 1649
extwar=sccs$v892,    # similar to 1650
himilexp=(sccs$v899==1)*1,
plunder=sccs$v912,
AP1=sccs$v921,           ###agricultural potential 1
AP2=sccs$v928,           ###agricultural potential 2
pathstress=sccs$v1260,
war=sccs$v1648,     # overall -- sum of internal and external
intwarB=sccs$v1649,  # similar to v891
extwarB=sccs$v1650,  # similar to v892
foodscarc=sccs$v1685,
sexratio=1+(sccs$v1689>85)+(sccs$v1689>115),
wagelabor=sccs$v1732,
CVrain=sccs$v1914/sccs$v1913   #no comma
)  
indep_vars<-c(
"famsize","exogamy",    "money","popdens","malesexag","ndrymonth","gath","hunt","fish",
"anim","brideprice","nuclearfam","ncmallow","cultints","tree","roots","cereals","settype","localjh",
"superjh","moralgods","segadlboys","plow","pigs","bovines","milk","agrlateboy","valchild","fratgrpstr",
"Whyte577","Whyte580","Whyte584","Whyte585","Whyte595","Whyte602","Whyte615","Whyte620","Whyte626","Whyte629",
"Whyte630","Whyte631","Whyte632","Whyte633","Whyte635","Paige657","femproduceND","Paige659","Paige660",
"Paige661","Paige662","fempower","interperviol","migr","Sanday664","Sanday665","Sanday666","Sanday667",
"Sanday668","Sanday669","Whyte718","Whyte719","Whyte720","Whyte721","Whyte722","Whyte723","Whyte724",
"Whyte725","Rohner798","Rohner799","Rohner800","Rohner801","Rohner802","Rohner803","Rohner804","Rohner805",
"Rohner806","Rohner807","Rohner808","Rohner809","Rohner810","Rohner811","Rohner812","Rohner813","foodtrade","fem_agri",
"dateobs","rain","temp","ecorich","pctFemPolyg","marrcaptives","femsubs","intwar","extwar","himilexp","plunder","AP1","AP2",
"pathstress","war","intwarB","extwarB","foodscarc","sexratio","wagelabor","CVrain"
)


restrict_vars=c("fratgrpstr","plow","fem_agri","marrcaptives","plunder","femproduceND","ecorich")
library(foreign)
#--Read in two weight matrices--
Wll<-as.matrix(read.dta("./examples/data/langwm.dta")[,-1])
Wdd<-as.matrix(read.dta("./examples/data/dist25wm.dta")[,c(-1,-2,-189)])
load("./examples/data/vaux.Rdata",.GlobalEnv)
my_aux = vaux
row.names(my_aux)<-NULL
#--remove the society name field--
my_aux<-my_aux[,-28]
name<-"general polygyny"
alias<-"DRWpolygyny2"
model=list(name=name,
           alias=alias,
           data=my_sccs,
           aux_data=my_aux,
           prox_list=list(language=Wll,distance=Wdd),
           dep_var="dep_var",
           indep_vars=indep_vars,
           restrict_vars=restrict_vars)
save(model,file=paste(alias,".Rdata",sep=""))
source("examples/src/run_model.R") #does for this model multiple imputation, two stage ols, saves to file to working directory.
lat<-sccs$v833.1
lon<-sccs$v833.2
 plot(lon,lat, cex=.1)
 ztxt=as.character(depvar) #above "my_sccs( " -- a new line inserted: depvar= as defined in my_sccs(  
 ztxt<-gsub("NaN",".",ztxt)
 text(lon,lat,ztxt)
ols_stats$restrict_stats
ols_stats$r2
ols_stats$restrict_diagnostics

1B Output CreateModelDRWpolygyny.R

                coef   Fstat       ddf pvalue    VIF
(Intercept)   3.2928  6.2110  193.5827 0.0135     NA
language     -0.7168  2.2763 9414.7259 0.1314 2.3153
distance      0.6272 13.3216 2037.9760 0.0003 2.1638
fratgrpstr    0.2503 12.9665   71.0921 0.0006 1.2979
plow         -0.5205  4.3369  229.4862 0.0384 1.4585
fem_agri      0.0085  4.7692  359.4181 0.0296 1.3843
marrcaptives  0.3179  9.3854   61.0798 0.0033 1.4399
plunder      -0.4536  6.2962  681.4361 0.0123 1.2654
femproduceND -0.4837  4.8211  126.4589 0.0299 1.0647
ecorich       0.3259  7.6741 1460.9834 0.0057 1.1819
>  ols_stats$r2
 R2:final model R2:IV_ language R2:IV_ distance 
      0.4637422       0.9964540       0.9979027 
>  ols_stats$restrict_diagnostics
               Fstat         df pvalue
RESET          0.111   3173.323  0.739
Wald.on.restrs 7.818      3.486  0.057
NCV            2.351   1926.180  0.125
SW.normal      0.119    871.410  0.730
lag..language  0.071 420589.855  0.790
lag..distance  0.080   2081.777  0.777

2A Copy 1A here and change your dep_var; get it running, then copy to ==3A Work on your indep_vars and restrict_vars==

setwd("C:/My Documents/sccs")

library(sccs)
data(sccs)                                 
depvar=sccs$v160
my_sccs<-data.frame(
dep_var=sccs$v160,
socname=sccs$socname,socID=sccs$"sccs#",
famsize=sccs$v68,    # similar to v80
exogamy=sccs$v72,
#famsizeB=sccs$v80,   # similar to v68
money=sccs$v155,
popdens=sccs$v156,
premarsexatt=sccs$v165,
premarsexfrq=sccs$v166,
malesexag=sccs$v175,
ndrymonth=sccs$v196,
gath=sccs$v203,
hunt=sccs$v204,
fish=sccs$v205,
anim=sccs$v206,
brideprice=(sccs$v208==1)*1,
nuclearfam=(sccs$v210<=3)*1,
ncmallow=sccs$v227,
cultints=sccs$v232,
tree=(sccs$v233==4)*1,
roots=(sccs$v233==5)*1,
cereals=(sccs$v233==6)*1,
settype=sccs$v234,
localjh=sccs$v236-1,
superjh=sccs$v237,
moralgods=sccs$v238,
segadlboys=sccs$v242,
plow=(sccs$v243>1)*1,
pigs=(sccs$v244==2)*1,
bovines=(sccs$v244==7)*1,
milk=(sccs$v245>1)*1,
agrlateboy=sccs$v300,
valchild=(sccs$v473+sccs$v474+sccs$v475+sccs$v476),
fratgrpstr=sccs$v570,
Whyte577=sccs$v577,    #take care using Whyte variables - only coded 1/2 the sample
Whyte580=sccs$v580,
Whyte584=sccs$v583,
Whyte585=sccs$v584,
Whyte595=sccs$v595,
Whyte602=sccs$v602,
Whyte615=sccs$v615,
Whyte620=sccs$v620,
Whyte626=sccs$v626,
Whyte629=sccs$v629,
Whyte630=sccs$v630,
Whyte631=sccs$v631,
Whyte632=sccs$v632,
Whyte633=sccs$v633,
Whyte635=sccs$v635,
#                    #take care using Paige variables - coded less than 1/2 the sample
Paige657=sccs$v657,  # summed in v663#  Paige657 Paige658 Paige659 Paige660 Paige661 Paige662
femproduceND=sccs$v658, #Paige658=sccs$v658,  # summed in v663 
Paige659=sccs$v659,  # summed in v663
Paige660=sccs$v660,  # summed in v663
Paige661=sccs$v661,  # summed in v663
Paige662=sccs$v662,  # summed in v663
fempower=sccs$v663, # sum of v657-662
interperviol=sccs$v666,   ###synonyms: violence / interviol / freintovio
migr=(sccs$v677==2)*1,
#                    #take care using Sanday variables - only coded 1/2 the sample
Sanday664=sccs$v664,  # summed in v669
Sanday665=sccs$v665,  # summed in v669
Sanday666=sccs$v666,  # summed in v669
Sanday667=sccs$v667,  # summed in v669
Sanday668=sccs$v668,  # summed in v669
Sanday669=sccs$v669, # sum of v664-668 
 #WHYTE Data Quality Whyte718 Whyte719 Whyte720 Whyte721 Whyte722 Whyte723 Whyte724 Whyte725
Whyte718=sccs$v718,    #take care using Whyte variables - only coded 1/2 the sample
Whyte719=sccs$v719,
Whyte720=sccs$v720,
Whyte721=sccs$v721,
Whyte722=sccs$v722,
Whyte723=sccs$v723,
Whyte724=sccs$v724,
Whyte725=sccs$v725,
 #Rohner Data Quality Codes 
Rohner798=sccs$v798,
Rohner799=sccs$v799,
Rohner800=sccs$v800,
Rohner801=sccs$v801,
Rohner802=sccs$v802,
Rohner803=sccs$v803,
Rohner804=sccs$v804,
Rohner805=sccs$v805,
Rohner806=sccs$v806,
Rohner807=sccs$v807,
Rohner808=sccs$v808,
Rohner809=sccs$v809,
Rohner810=sccs$v810,
Rohner811=sccs$v811,
Rohner812=sccs$v812,
Rohner813=sccs$v813,
foodtrade=sccs$v819,
fem_agri=sccs$v821, 
dateobs=sccs$v838,
rain=sccs$v854,
temp=sccs$v855,
#ecorich=sccs$v857,
ecorich=(sccs$v857==3|sccs$v857==4)*1+(sccs$v857==5)*2,
pctFemPolyg=sccs$v872,
marrcaptives=sccs$v870,
femsubs=sccs$v890,
intwar=sccs$v891,    # similar to 1649
extwar=sccs$v892,    # similar to 1650
himilexp=(sccs$v899==1)*1,
plunder=sccs$v912,
AP1=sccs$v921,           ###agricultural potential 1
AP2=sccs$v928,           ###agricultural potential 2
pathstress=sccs$v1260,
war=sccs$v1648,     # overall -- sum of internal and external
intwarB=sccs$v1649,  # similar to v891
extwarB=sccs$v1650,  # similar to v892
foodscarc=sccs$v1685,
sexratio=1+(sccs$v1689>85)+(sccs$v1689>115),
wagelabor=sccs$v1732,
CVrain=sccs$v1914/sccs$v1913   #no comma
)  
indep_vars<-c(
"famsize","exogamy",    "money","popdens","malesexag","ndrymonth","gath","hunt","fish",
"anim","brideprice","nuclearfam","ncmallow","cultints","tree","roots","cereals","settype","localjh",
"superjh","moralgods","segadlboys","plow","pigs","bovines","milk","agrlateboy","valchild","fratgrpstr",
"Whyte577","Whyte580","Whyte584","Whyte585","Whyte595","Whyte602","Whyte615","Whyte620","Whyte626","Whyte629",
"Whyte630","Whyte631","Whyte632","Whyte633","Whyte635","Paige657","femproduceND","Paige659","Paige660",
"Paige661","Paige662","fempower","interperviol","migr","Sanday664","Sanday665","Sanday666","Sanday667",
"Sanday668","Sanday669","Whyte718","Whyte719","Whyte720","Whyte721","Whyte722","Whyte723","Whyte724",
"Whyte725","Rohner798","Rohner799","Rohner800","Rohner801","Rohner802","Rohner803","Rohner804","Rohner805",
"Rohner806","Rohner807","Rohner808","Rohner809","Rohner810","Rohner811","Rohner812","Rohner813","foodtrade","fem_agri",
"dateobs","rain","temp","ecorich","pctFemPolyg","marrcaptives","femsubs","intwar","extwar","himilexp","AP1","AP2",
"pathstress","war","intwarB","extwarB","foodscarc","sexratio","wagelabor","CVrain"
)
restrict_vars=c("famsize","exogamy","money","moralgods","nuclearfam","sexratio","wagelabor",)
library(foreign)
#--Read in two weight matrices--
Wll<-as.matrix(read.dta("./examples/data/langwm.dta")[,-1])
Wdd<-as.matrix(read.dta("./examples/data/dist25wm.dta")[,c(-1,-2,-189)])
load("./examples/data/vaux.Rdata",.GlobalEnv)
my_aux = vaux
row.names(my_aux)<-NULL
#--remove the society name field--
my_aux<-my_aux[,-28]
name<-"general polygyny"
alias<-"DRWpolygyny2"
model=list(name=name,
           alias=alias,
           data=my_sccs,
           aux_data=my_aux,
           prox_list=list(language=Wll,distance=Wdd),
           dep_var="dep_var",
           indep_vars=indep_vars,
           restrict_vars=restrict_vars)
save(model,file=paste(alias,".Rdata",sep=""))
source("examples/src/run_model.R") #does for this model multiple imputation, two stage ols, saves to file to working directory.
lat<-sccs$v833.1
lon<-sccs$v833.2
plot(lon,lat, cex=.1)
ztxt=as.character(depvar) #above "my_sccs( " -- a new line inserted: depvar= as defined in my_sccs( 
text(lon,lat,ztxt)
ols_stats$restrict_stats
ols_stats$r2
ols_stats$restrict_diagnostics
table(depvar)

2B Your 2A Results: Give you dep_var name and vNumber in sccs

             coef Fstat        ddf pvalue   VIF

(Intercept) 0.230 0.199 5154.036 0.656 NA language -0.066 0.009 1586.065 0.926 2.179 distance 0.875 6.702 1312.956 0.010 2.108 cultints -0.046 0.674 135553.366 0.412 1.580 bovines 0.366 3.062 582154.648 0.080 1.575 ncmallow -0.017 0.261 3944.007 0.610 1.081 gath -0.032 0.254 11365.810 0.615 1.526 malesexag 0.026 0.213 33.714 0.648 1.239 > ols_stats$r2

R2:final model R2:IV_ language R2:IV_ distance 
    0.09249096      0.99478393      0.99530588 

> ols_stats$restrict_diagnostics

               Fstat         df pvalue

RESET 0.417 391.380 0.519 Wald.on.restrs 3.391 8.955 0.099 NCV 0.732 530.184 0.393 SW.normal 27.637 14433.564 0.000 lag..language 0.016 77320.831 0.900 lag..distance 0.068 266892.210 0.795 > table(depvar) depvar

1  2  3  4 

11 40 6 9 >

3A Work on your indep_vars and restrict_vars