-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathBrassRelationalLogitCode.R
113 lines (95 loc) · 5.67 KB
/
BrassRelationalLogitCode.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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
##############################################################################################################################
##R CODE FOR THE BRASS RELATIONAL LOGIT MODEL OF MORTALITY
##
##EDDIE HUNSINGER, FEBRUARY 2011 (LAST UPDATED OCTOBER 2022)
##http://www.demog.berkeley.edu/~eddieh/
##edyhsgr@gmail.com
##
##EXAMPLE DATA IS LINKED, SO YOU SHOULD BE ABLE TO SIMPLY COPY ALL AND PASTE INTO R
##
##A JUPYTER NOTEBOOK (PYTHON) VERSION OF THE CODE IS AVAILABLE AT:
## https://github.com/AppliedDemogToolbox/Hunsinger_BrassRelationalLogit/blob/master/BrassRelationalLogitCode.ipynb
##
##THERE IS NO WARRANTY FOR THIS CODE
##
##A PYTHON NOTEBOOK VERSION OF THIS PROCESS IS AVAILABLE AT https://github.com/AppliedDemogToolbox/Hunsinger_BrassRelationalLogit/blob/master/BrassRelationalLogitCode.ipynb
##############################################################################################################################
##############################################################################################################################
#STEP 1: Read in and plot the survivorship data (lxSample is collected data, lxSSA is Social Security Administration data, to
#be used as the standard-- so the goal is to adjust (the more reliable) lxSSA to meet the broad shape and level of lxSample)
##############################################################################################################################
lxSample<-read.table(file="https://raw.githubusercontent.com/AppliedDemogToolbox/Hunsinger_BrassRelationalLogit/master/lxSample.csv",header=TRUE,sep=",")
lxSSA<-read.table(file="https://raw.githubusercontent.com/AppliedDemogToolbox/Hunsinger_BrassRelationalLogit/master/lxSSA.csv",header=TRUE,sep=",")
plot(lxSample$X2000f,type="l",
main="Brass relational logit model demonstration",
xlab="age",
ylab="survivorship",
col="blue",lwd=5,axes=FALSE)
axis(side=2,las=2)
axis(side=1,at=1:length(lxSample$x),
labels=lxSample$x,las=2)
points(lxSSA$X2000f,type="l",lty=2,col="red",lwd=5)
legend(x="left",legend=c("sample","standard"),
col=c("blue","red"),lty=c(1,2),lwd=c(5,5),bg="white")
##############################################################################################################################
#STEP 2: Enter the Brass Relational Logit function
##############################################################################################################################
Brass<-function(lx,lxBase){
Yx<-.5*log(lx/(1-lx))
YxBase<-.5*log(lxBase/(1-lxBase))
Estimate<-lm(Yx[2:length(lxBase)]~YxBase[2:length(lxBase)])
Alpha<-Estimate$coefficients[1]
Beta<-Estimate$coefficients[2]
Brass<- data.frame(Beta=Beta,Alpha=Alpha)
return(Brass)}
##############################################################################################################################
#STEP 3: Using the Brass function, estimate the coefficients for adjusting lxSSA
##############################################################################################################################
BrassCoefficients<-Brass(lxSample$X2000f,lxSSA$X2000f)
BrassCoefficients
##############################################################################################################################
#STEP 4: Adjust lxSSA with the estimated coefficients and plot the adjusted lxSSA
##############################################################################################################################
Yx<-NULL
for (i in 1:length(lxSSA$X2000f)){Yx[i]<-.5*log(lxSSA$X2000f[i]/(1-lxSSA$X2000f[i]))}
lxSSAAdjusted<-NULL
for (i in 1:length(lxSSA$X2000f)){lxSSAAdjusted[i]<-1/(1+exp(-2*BrassCoefficients$Alpha-2*BrassCoefficients$Beta*Yx[i]))}
points(lxSSAAdjusted,type="l",lty=3,col="purple",lwd=5)
legend(x="left",legend=c("sample","standard","adjusted standard"),
col=c("blue","red","purple"),lty=c(1,2,3),lwd=c(5,5,5),bg="white")
##############################################################################################################################
#STEP 5: Plot some effects of change in Brass Alpha and Brass Beta
##############################################################################################################################
lxSSAHighAlpha<-NULL
for (i in 1:length(lxSSA$X2000f)){lxSSAHighAlpha[i]<-1/(1+exp(-2*(.5)-2*(1)*Yx[i]))}
lxSSALowAlpha<-NULL
for (i in 1:length(lxSSA$X2000f)){lxSSALowAlpha[i]<-1/(1+exp(-2*(-.5)-2*(1)*Yx[i]))}
lxSSAHighBeta<-NULL
for (i in 1:length(lxSSA$X2000f)){lxSSAHighBeta[i]<-1/(1+exp(-2*(0)-2*(1.5)*Yx[i]))}
lxSSALowBeta<-NULL
for (i in 1:length(lxSSA$X2000f)){lxSSALowBeta[i]<-1/(1+exp(-2*(0)-2*(.5)*Yx[i]))}
plot(lxSSA$X2000f,type="l",
main="Brass relational logit model demonstration",
xlab="age",
ylab="survivorship",
col="purple",lwd=5,axes=FALSE)
axis(side=2,las=2)
axis(side=1,at=1:length(lxSample$x),
labels=lxSample$x,las=2)
points(lxSSAHighAlpha,type="l",lty=2,col="red",lwd=5)
points(lxSSALowAlpha,type="l",lty=3,col="blue",lwd=5)
legend(x="left",legend=c("standard","standard with 0.5 alpha","standard with -0.5 alpha"),
col=c("purple","red","blue"),lty=c(1,2,3),lwd=c(5,5,5),bg="white")
plot(lxSSA$X2000f,type="l",
main="Brass relational logit model demonstration",
xlab="age",
ylab="survivorship",
col="purple",lwd=5,axes=FALSE)
axis(side=2,las=2)
axis(side=1,at=1:length(lxSample$x),
labels=lxSample$x,las=2)
points(lxSSAHighBeta,type="l",lty=2,col="red",lwd=5)
points(lxSSALowBeta,type="l",lty=3,col="blue",lwd=5)
legend(x="left",legend=c("standard","standard with 1.5 beta","standard with 0.5 beta"),
col=c("purple","red","blue"),lty=c(1,2,3),lwd=c(5,5,5),bg="white")
#write.table(###, file="G:/###/###.csv", sep=",")