-
Notifications
You must be signed in to change notification settings - Fork 0
/
regress.R
72 lines (71 loc) · 4.66 KB
/
regress.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
decimal1_2<-function(x){ # for conversion of character variables to two decimal places.
if (nchar(x)==1) paste(x,".00",sep="")
else if (nchar(x)==3) paste(x,"0",sep="")
else return(x)
}
df2<-subset(df,previoustb!="Missing")
df2$obs<-NULL
df2$silicosis2<-0
df2$silicosis2[df2$silicosis=="probable"]<-1
m1<-glm(tb2~silicosis2*country2+factor(agecat3)+dwelling+grade+previoustb+mine,data=df2,family="binomial")
orm1<-or_wrapper(m1)
# with interactions taken into account
df2$agecat3<-factor(df2$agecat3)
#print("here?")
d<-datadist(df2)
#print("before this?")
options(datadist='d')
a<-lrm(tb2~silicosis2*country2+agecat3+dwelling+grade+previoustb+mine,data=df2)
#a
#summary(a)
sil_L<-as.character(round(summary(a,country2="Lesotho")[2,4],2))
sil_M<-as.character(round(summary(a,country2="Mozambique")[2,4],2))
sil_SA<-as.character(round(summary(a,country2="South Africa")[2,4],2))
sil_L_lci<-as.character(round(summary(a,country2="Lesotho")[2,6],2))
sil_M_lci<-as.character(round(summary(a,country2="Mozambique")[2,6],2))
sil_SA_lci<-as.character(round(summary(a,country2="South Africa")[2,6],2))
sil_L_uci<-as.character(round(summary(a,country2="Lesotho")[2,7],2))
sil_M_uci<-as.character(round(summary(a,country2="Mozambique")[2,7],2))
sil_SA_uci<-as.character(round(summary(a,country2="South Africa")[2,7],2))
sil_L_p<-paste(decimal1_2(sil_L)," (",decimal1_2(sil_L_lci),"-",decimal1_2(sil_L_uci),")",sep="")
sil_M_p<-paste(decimal1_2(sil_M)," (",decimal1_2(sil_M_lci),"-",decimal1_2(sil_M_uci),")",sep="")
sil_SA_p<-paste(decimal1_2(sil_SA)," (",decimal1_2(sil_SA_lci),"-",decimal1_2(sil_SA_uci),")",sep="")
nice_model<-data.frame("R"=as.character(),
"I"=as.character(),
"O"=as.character())
nice_model<-rbind(nice_model,c( paste(" n = ",as.character(length(df2$idnum)),sep=""),"","" ))
names(nice_model)<-c("R","I","O")
nice_model$R<-as.character(nice_model$R)
nice_model$I<-as.character(nice_model$I)
nice_model$O<-as.character(nice_model$O)
nice_model<-rbind(nice_model,c("Silicosis","Lesotho",sil_L_p))
nice_model$R<-as.character(nice_model$R)
nice_model$I<-as.character(nice_model$I)
nice_model$O<-as.character(nice_model$O)
nice_model<-rbind(nice_model,c("","Mozambique",sil_M_p))
nice_model<-rbind(nice_model,c("","South Africa",sil_SA_p))
nice_model<-rbind(nice_model,c("Age <35","",""))
nice_model<-rbind(nice_model,c("Age 35-44","",paste(round(orm1[5,1],2)," (",round(orm1[5,2],2),"-",round(orm1[5,3],2),")",sep="")))
nice_model<-rbind(nice_model,c("Age 45+","",paste(round(orm1[6,1],2)," (",round(orm1[6,2],2),"-",round(orm1[6,3],2),".00",")",sep="")))
nice_model<-rbind(nice_model,c("Accomodation hostel","",""))
nice_model<-rbind(nice_model,c("Accomodation non-hostel","",paste(round(1/(orm1[7,1]),2)," (",round(1/(orm1[7,3]),2),"-",round(1/(orm1[7,2]),2),")",sep=""))) #Need to set this to 1/x due to levels of factor.
nice_model<-rbind(nice_model,c("Job grade skilled","",""))
nice_model<-rbind(nice_model,c("Job grade unskilled","",paste(round(orm1[8,1],2)," (",round(orm1[8,2],2),"-",round(orm1[8,3],2),")",sep="")))
nice_model<-rbind(nice_model,c("No previous TB","",""))
nice_model<-rbind(nice_model,c("Previous TB","",paste(round(orm1[9,1],2)," (",round(orm1[9,2],2),"-",round(orm1[9,3],2),")",sep="")))
nice_model<-rbind(nice_model,c("Mine A","",""))
nice_model<-rbind(nice_model,c("Mine B","",paste(round(orm1[10,1],2)," (",round(orm1[10,2],2),"-",round(orm1[10,3],2),")",sep="")))
#nice_model<-rbind(nice_model,c("South Africa","No/possible silicosis",""))
#nice_model<-rbind(nice_model,c("","Probable silicosis",""))
L_s0<-paste(round(summary(a,silicosis2=0)[6,4],2)," (",round(summary(a,silicosis2=0)[6,6],2),"-",round(summary(a,silicosis2=0)[6,7],2),")",sep="")
M_s0<-paste(round(summary(a,silicosis2=0)[4,4],2)," (",round(summary(a,silicosis2=0)[4,6],2),"-",round(summary(a,silicosis2=0)[4,7],2),")",sep="")
L_s1<-paste(round(summary(a,silicosis2=1)[6,4],2)," (",round(summary(a,silicosis2=1)[6,6],2),"-",round(summary(a,silicosis2=1)[6,7],2),")",sep="")
M_s1<-paste(round(summary(a,silicosis2=1)[4,4],2)," (",round(summary(a,silicosis2=1)[4,6],2),"-",round(summary(a,silicosis2=1)[4,7],2),")",sep="")
#nice_model<-rbind(nice_model,c("Lesotho","No/possible silicosis",L_s0))
#nice_model<-rbind(nice_model,c("","Probable silicosis",L_s1))
#nice_model<-rbind(nice_model,c("Mozambique","No/possible silicosis",M_s0))
#nice_model<-rbind(nice_model,c("","Probable silicosis",M_s1))
#summary(a,silicosis2=0)
#summary(a,silicosis2=1)
names(nice_model)<-c("Risk factor","Interacting level","Odds ratio (95% CI)")
reg_out<-xtable(nice_model,caption="Multivariable logistic regression analysis of silicosis as a risk factor for tuberculosis in male mine workers in South Africa, 2006.",label="reg")