-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathDewarManual.R
63 lines (48 loc) · 3.46 KB
/
DewarManual.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#################################################################
#
# R-Studio @ Macquarie University, Sydney, Australia
# 15.05.2017 - Adela Sobotkova
#
# Topic: Dealing with Contemporaneity in Survey data through a Census approach (a la Dewar 1991)
#
################################################################
#### Set Workspace, if your data is tabular
#setwd("D:/Users/MQ20149304/Documents/RStudio/Workshop1")
# to be continued >> Work in Progress, ask Adela about DewarUpgrade
###If no workspace needed, enter parameters directly here (currently Shepelah region in Nitsan's dataset)
phase=207 # length in years of $ phase
vill=32 #product of a+b in Dewar
settle=39# (c+d)/phaselength, site establishment rate * 1000 (to match random integer data generated by sample function)
aban=58 # (a+d)/phaselength,site abandonment rate * 1000 (to match random integer data generated by sample function)
year=1 #starting year of the phase iteration
### Create a matrix in which you will aggregate the interim data
phase.results<-matrix(year, nrow = 0, ncol = 3,byrow=T)
dimnames(phase.results)<-list(
c(), # row names
c("Year", "CurrentVillages", "SD")) # column names
print(phase.results) # just a test message that this part of code is not broken, it shoud only print the column names.
### Create a loop that models the emergence of new villages on the basis of the existing site abandonment and creation rates
## The loop and random number generation
while(year<=phase){iter <- sample(1:1000,999,rep=T) # creates 1000 random numbers for every year of the phase, runs 999 times
iter2 <- ifelse(iter<=settle,1,0) # checks if the random number is below the settlement rate and if so, output 1, signalling that a settlement was newly inhabited
iter3 <- ifelse(iter<=aban,1,0) # checks if the random number is below the abandonment rate and if so, output 1, signalling that a village was abandoned
new<-iter2+vill-iter3 # adds the newly settled and substracts the abandoned settlements
vill<-mean(new) # calculates the mean of the new villages in a given year on the basis of 999 iterations
x<-sd(new) # calculates the standard deviation of new villages in a given year
year.results<-c(year,vill,x) #creates a vector of data from the annual statistic
phase.results<-rbind(phase.results,c(year,vill,x)) #combines all the individual annual iter into a table for the phase
year<-year+1
}
### Results #####
Occ<-mean(phase.results[,2]) # number of simultaneously occupied villages
Occ.sd<-sd(phase.results[,2]) # the standard deviation in the simultaneously occupied villages
Use.span<-1/(aban*0.001/Occ) # mean use span of the villages in this phase. the multiplication reconverts the abandonment rate to original real number
Low.span<-1/(aban*0.001/(Occ-Occ.sd)) # mean use span minus 1 sd. the multiplication reconverts the abandonment rate to original real number
Hi.span<-1/(aban*0.001/(Occ+Occ.sd)) # mean use span plus 1 sd. the multiplication reconverts the abandonment rate to original real number
Span.err<- (Hi.span-Low.span)/2 # the standard deviation in the mean use span
### Final Messages #####
message("During this phase the number of simultaneously occupied villages is ca ",Occ,"+/-",Occ.sd)
message("The mean use span with 1 standard deviation is ",Use.span,"+/-",Span.err)
message("Less one sd the span estimate is ",Low.span)
message("Plus one sd the span estimate is ",Hi.span)
##### End of the Exercise ####################