FinalProject
Anthony Colombo
11/7/2020
Library
library(magrittr)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(reshape2)
library(magrittr);library(dplyr)
library(lattice)Data management
we add informative labels to the factors. A basic EDA showed that age and doctor education did not vary across misuse, but temperature, self medication, and wrong diagnosis did vary.
anti<-read.csv("~/Documents/school/PM511c/final-project/antibiotics.csv")
og<-anti
anti$dich_misuse<-ifelse(anti$misuse==1,0,1)
anti$misuse<-as.character(anti$misuse)
anti$misuse[which(anti$misuse=="1")]<-'correct use'
anti$misuse[which(anti$misuse=="2")]<-'misuse of 1 antibiotic'
anti$misuse[which(anti$misuse=="3")]<-'misuse of several antibiotics'
anti$misuse<-factor(anti$misuse)
## adding labels
anti$Paymed<-factor(ifelse(anti$Paymed==1,"Patient self paying","Patient not self paying"))
anti$Selfmed<-factor(ifelse(anti$Selfmed==1,"Patient self medicating","Patient not self medicating"))
anti$Wrdiag<-factor(ifelse(anti$Wrdiag==1,"Wrong diagnosis","Correct diagnosis"))
anti$DRed<-factor(anti$DRed)
anti$WHO<-factor(ifelse(anti$WHO==1,"WHO participation","WHO non-participation"))
## EDA
# aggregate(age~misuse,og,FUN=mean)
# aggregate(temp~misuse,og,FUN=mean)
# aggregate(Selfmed~misuse,og,FUN=mean)
# aggregate(Wrdiag~misuse,og,FUN=mean)
# aggregate(DRed~misuse,og,FUN=mean)
# aggregate(WHO~misuse,og,FUN=mean)Table 1
The table 1 is created using the annotated labels across each factor.
library(table1)##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<- anti.table<-anti
colnames(anti.table)<-c("Hospital",
"Doctor.Identifier",
"Patient_Self_Payment",
"Age",
"Temperature",
"Prior_self_medication",
"Antibiotic_misuse",
"Wrong_diagnosis",
"Doctor_education_Level",
"WHO_participation",
"misuse_category")
anti.table[,'misuse_category']<-ifelse(anti.table[,'misuse_category']==0,"Correct antibiotic use","Misuse of at least 1 antibiotic")
anti.table[,'misuse_category']<-factor( anti.table[,'misuse_category'])
T1<-table1(~Patient_Self_Payment+Age+Temperature+Antibiotic_misuse+Prior_self_medication+Wrong_diagnosis+Doctor_education_Level+WHO_participation|misuse_category, anti.table)Cohort summary measures of level 2 and 3
- The proportion of hospitals that participated in WHO is 27.8%. we count the number of children per hospital.
## Hospital summary
dd<-(table(anti.table$Hospital,anti.table$WHO_participation))
print(head(dd))##
## WHO non-participation WHO participation
## 1 62 0
## 2 15 0
## 3 46 0
## 4 45 0
## 5 31 0
## 6 30 0 dd[,1]<-ifelse(dd[,1]!=0,1,0)
dd[,2]<-ifelse(dd[,2]!=0,1,0)
100*(10/36)## [1] 27.77778Proportion of children per doctor (level 2 summary)
- The proportion of self-taught doctors, among doctor relative proportions is 1.5%. The medical doctor (6) relative proportion is 9.7%.
## Doctor level summary
dd<-(table(anti.table$Doctor.Identifier,anti.table$Doctor_education_Level))
dd<-(apply(dd,2,function(x) ifelse(x!=0,1,0)))
100*colSums(dd)/sum(dd)## 1 2 3 4 5 6
## 1.492537 29.104478 8.955224 44.776119 5.970149 9.701493Simple linear regression
- After adjusting for WHO participation, Age is not different among usage (p=0.0502). Temperature does differ, further temperature is not associated with prior self medication (p=0.94).
##age comparison differences
summary(lm(Age~WHO_participation+misuse_category,anti.table))##
## Call:
## lm(formula = Age ~ WHO_participation + misuse_category, data = anti.table)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3584 -1.1623 -0.1623 1.1316 2.3277
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 2.16230 0.08794 24.588
## WHO_participationWHO participation -0.49000 0.13038 -3.758
## misuse_categoryMisuse of at least 1 antibiotic 0.19610 0.10003 1.960
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## WHO_participationWHO participation 0.000183 ***
## misuse_categoryMisuse of at least 1 antibiotic 0.050265 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.324 on 852 degrees of freedom
## Multiple R-squared: 0.02455, Adjusted R-squared: 0.02226
## F-statistic: 10.72 on 2 and 852 DF, p-value: 2.518e-05 confint(lm(Age~WHO_participation+misuse_category,anti.table))## 2.5 % 97.5 %
## (Intercept) 1.9896926983 2.3349029
## WHO_participationWHO participation -0.7459034683 -0.2340881
## misuse_categoryMisuse of at least 1 antibiotic -0.0002273931 0.3924322 ### Temperature comparison
summary(lm(Temperature~WHO_participation+misuse_category,anti.table))##
## Call:
## lm(formula = Temperature ~ WHO_participation + misuse_category,
## data = anti.table)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.07624 -0.68285 0.01715 0.72314 2.41716
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.98285 0.06272 31.617
## WHO_participationWHO participation 0.09339 0.09298 1.004
## misuse_categoryMisuse of at least 1 antibiotic -0.51197 0.07134 -7.177
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## WHO_participationWHO participation 0.315
## misuse_categoryMisuse of at least 1 antibiotic 1.55e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9439 on 852 degrees of freedom
## Multiple R-squared: 0.06288, Adjusted R-squared: 0.06068
## F-statistic: 28.58 on 2 and 852 DF, p-value: 9.662e-13 ##
summary(lm(Temperature~WHO_participation+Prior_self_medication,anti.table))##
## Call:
## lm(formula = Temperature ~ WHO_participation + Prior_self_medication,
## data = anti.table)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8275 -0.8106 -0.1049 0.7894 2.7951
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.610555 0.041327 38.971
## WHO_participationWHO participation 0.216978 0.094224 2.303
## Prior_self_medicationPatient self medicating -0.005666 0.076152 -0.074
## Pr(>|t|)
## (Intercept) <2e-16 ***
## WHO_participationWHO participation 0.0215 *
## Prior_self_medicationPatient self medicating 0.9407
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.972 on 852 degrees of freedom
## Multiple R-squared: 0.006232, Adjusted R-squared: 0.003899
## F-statistic: 2.671 on 2 and 852 DF, p-value: 0.06974Proportion of WHO comparision
- There is 17% less misuse across hospitals which participated in WHO compared to those that did not.
x<-as.data.frame.matrix(table(anti.table$Hospital,anti.table$misuse_category))
x$ID<-rownames(x)
x$ID<-factor(x$ID,levels=rownames(x))
x$WHO<-anti.table[match(x$ID,anti.table$Hospital),'WHO_participation']
### Proportion of WHO
x$prop<-x[,2]/(x[,1]+x[,2])
summary(lm(prop~WHO,x))##
## Call:
## lm(formula = prop ~ WHO, data = x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.299006 -0.066322 0.001587 0.089877 0.236708
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.72074 0.02756 26.152 < 2e-16 ***
## WHOWHO participation -0.17174 0.05229 -3.284 0.00237 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1405 on 34 degrees of freedom
## Multiple R-squared: 0.2408, Adjusted R-squared: 0.2185
## F-statistic: 10.79 on 1 and 34 DF, p-value: 0.002375Bullet 1b
- we compute the variability across doctor and hospital by looking at the residual ICC in each variance compoenent model.
- the ICC of doctor : 0.15
- The ICC of hospital : 0.104
- we manually compute the total variance across the misuse (binary factor) using the law of total variance \(V(Doctor)=E(Var(Doctor|I_{Misuse}))+V(E(Doctor|I_{Misuse}))\)
- The total variation of the number of pediatric patients visiting each doctor across correct/incorrect antibiotic use is 6.49.
- Similarly, the total variation of the number of patients visiting each hospital across correct/incorrect use is 108.35.
library(lme4);library(magrittr);library(dplyr)## Loading required package: Matrixlibrary(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(ggplot2)
library(hrbrthemes)## NOTE: Either Arial Narrow or Roboto Condensed fonts are required to use these themes.## Please use hrbrthemes::import_roboto_condensed() to install Roboto Condensed and## if Arial Narrow is not on your system, please see https://bit.ly/arialnarrow fit1 <- glmer(dich_misuse ~(1|doc), data = anti, family = binomial,
nAGQ=25,control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
ICC.doc<-(0.762)^2/(0.762^2+pi^2/3)
doc<-anti%>%group_by(factor(doc),dich_misuse)%>%summarise(N=n())%>%data.frame## `summarise()` regrouping output by 'factor(doc)' (override with `.groups` argument) fit2 <- glmer(dich_misuse ~(1|hosp), data = anti, family = binomial,
nAGQ=25,control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
ICC.hos<-(0.6188)^2/(0.6188^2+pi^2/3)
## add a graph here.
x<-as.data.frame.matrix(table(anti.table$Doctor.Identifier,anti.table$misuse_category))
x$ID<-rownames(x)
x$ID<-factor(x$ID,levels=rownames(x))
ggplot(melt(x),aes(x=ID,y=value,color=variable))+geom_point()+ylab("Number of antibiotic use/miuse")+scale_color_manual("Antibiotic use",values=c('blue','red'))+xlab("Doctor ID")+theme_ipsum(base_family='Arial')+theme(axis.text.x=element_text(angle=45,hjust=1,size=3))+ggtitle("Variability of misuse across doctors")## Using ID as id variables
## total variance across doctors of counts of patient misuse for 2 categories (correct or at least 1).
var_between=var(apply(x[,-3],2,function(x) mean(x)))
var_within=mean(apply(x[,-3],2,function(x) var(x)))
total_var_dr=var_between+var_within
## HOSPTIAL VARIATION
x<-as.data.frame.matrix(table(anti.table$Hospital,anti.table$misuse_category))
x$ID<-rownames(x)
x$ID<-factor(x$ID,levels=rownames(x))
ggplot(melt(x),aes(x=ID,y=value,color=variable))+geom_point()+ylab("Number of antibiotic use/miuse")+scale_color_manual("Antibiotic use",values=c('blue','red'))+xlab("Hospital ID")+theme_ipsum(base_family='Arial')+theme(axis.text.x=element_text(angle=45,hjust=1,size=5))+ggtitle("Variability of misuse across hospitals")## Using ID as id variables
# barplot(table(table(anti.table$Doctor.Identifier)),
# main='Distribution of numbers of pediatric visits across doctor',
# xlab='Number of doctor visits per doctor')
## total variance across hospitals of counts of patient misuse for 2 categories (correct or at least 1).
var_between=var(apply(x[,-3],2,function(x) mean(x)))
var_within=mean(apply(x[,-3],2,function(x) var(x)))
total_var=var_between+var_within
# barplot(table(table(anti.table$Hospital)),
# main='Distribution of numbers of visits across hospital',
# xlab='Number of doctor visits per hospital')Count summary
- We do observe major differences of the number of children visiting hospitals with increased misuse from hospitals 1,2,4, and 22
library(ComplexHeatmap)## Loading required package: grid## ========================================
## ComplexHeatmap version 2.4.3
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite:
## Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ======================================== ## HOSPTIAL VARIATION
x<-as.data.frame.matrix(table(anti.table$Hospital,anti.table$misuse_category))
x$ID<-rownames(x)
x$ID<-factor(x$ID,levels=rownames(x))
x$WHO<-anti.table[match(x$ID,anti.table$Hospital),'WHO_participation']
ha = rowAnnotation(WHO = x$WHO,
gp=gpar(col=ifelse(x$WHO=='WHO participation',"black","pink")))
heats<-Heatmap(x[,c(1,2)],name='Antibiotic counts',column_title='Antibiotic usage count per hospital',left_annotation=ha)## Warning: The input is a data frame, convert it to the matrix. draw(heats)
x$prop<-x[,2]/(x[,1]+x[,2])
summary(lm(prop~WHO,x))##
## Call:
## lm(formula = prop ~ WHO, data = x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.299006 -0.066322 0.001587 0.089877 0.236708
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.72074 0.02756 26.152 < 2e-16 ***
## WHOWHO participation -0.17174 0.05229 -3.284 0.00237 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1405 on 34 degrees of freedom
## Multiple R-squared: 0.2408, Adjusted R-squared: 0.2185
## F-statistic: 10.79 on 1 and 34 DF, p-value: 0.002375logistic regression unadjusted models
- Fitting a 3-level un-adjusted model for the association between antibiotic misuse and WHO participation, the AIC was 1022.6, with the WHO participation as significant (p=2.89e-03).
- fitting a 2-level model for just the hospital, had a AIC 1022.7.
- the 3-level model does not differ from only including the hospital as random effect (un-nested) (p=0.14, ANOVA), however the 3 level model had slightly lower AIC (1022.613) compared to only hospital effects (1022.712).
- The model is defined as i,j,k represents pediatric patient, doctors, and hospital. and \(Y_{ik}\) is an indicator for 1 for misuse of at least one antibiotic, or 0 correct use.
- \(logit(Pr(Y_{ik}==1))= (\beta_{0}+U_{0k}+U_{0jk})+\beta_{1}WHO_k\)
- The WHO 3-level model estimate yielded an estimate -1.021 (-1.58,-0.44).
## unadjusted WHO
fit3 <- glmer(misuse_category ~ WHO_participation+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
fit4 <- glmer(misuse_category ~ WHO_participation+(1|Hospital),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
anova(fit4,fit3) ## requires variable selection.
fit3.ci<-confint(fit3)## Computing profile confidence intervals ...Unadjusted WHO Odds interpretation.
- For a given doctor, in a given hospital, the unadjusted odds for antibiotic misuse of participants with at least 1 antibiotic is 0.36 (0.21,0.64) times lower in hospitals which participated in WHO program compared to hospitals which did not participate (p=2.9e-03).
exp(fit3.ci)## 2.5 % 97.5 %
## .sig01 1.000000 2.164086
## .sig02 1.000000 1.991890
## (Intercept) 2.239918 3.793625
## WHO_participationWHO participation 0.206143 0.643271 exp(-1.021)## [1] 0.3602345Adjusting for continuous features and selecting features in a model.
- The raw group averages of misuse across age, are similar using a naive comparison.
- Age was not significant as a fixed effect (p=0.11), nor interactions between child age and WHO (p=0.84). There was not a significant quadratic trend either (p=0.82).
- Because age was not a significant variable, we did not group center.
- We did grand centering on age (2.227) on age.
- Fitting temperature, we identified a quadratic polynomial trend for temperature (p=0.02), however the cubic term was not significant compared to only a quadratic feature for temperature (p=0.91). The interaction between WHO and temperature was marginally significant (p=0.07) which warrants future consideration, despite the main effect of temperature being not significant (p=0.88)!
- The AIC of the temperature model that included quadratic term is 966.3
aggregate(age~misuse,og,FUN=mean) ## do we center age?
fit.age <- glmer(misuse_category ~ WHO_participation*Age+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
fit.age.poly <- glmer(misuse_category ~ WHO_participation*Age+I(Age**2)+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
anova(fit.age,fit.age.poly)summary(fit.age.poly)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## misuse_category ~ WHO_participation * Age + I(Age^2) + (1 | Hospital/Doctor.Identifier)
## Data: anti.table
## Control:
## glmerControl(optimizer = "bobyqa", boundary.tol = 3e-05, optCtrl = list(maxfun = 3e+05))
##
## AIC BIC logLik deviance df.resid
## 1025.1 1058.4 -505.5 1011.1 848
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2515 -0.9663 0.5157 0.6300 1.2389
##
## Random effects:
## Groups Name Variance Std.Dev.
## Doctor.Identifier:Hospital (Intercept) 0.2033 0.4509
## Hospital (Intercept) 0.1395 0.3734
## Number of obs: 855, groups: Doctor.Identifier:Hospital, 134; Hospital, 36
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.84267 0.23527 3.582 0.000341 ***
## WHO_participationWHO participation -1.03780 0.42212 -2.459 0.013950 *
## Age 0.06073 0.21898 0.277 0.781536
## I(Age^2) 0.01149 0.04975 0.231 0.817409
## WHO_participationWHO participation:Age 0.03830 0.17000 0.225 0.821765
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WHO_pWHOp Age I(A^2)
## WHO_prtWHOp -0.408
## Age -0.716 0.136
## I(Age^2) 0.546 -0.031 -0.952
## WHO_pWHOp:A 0.302 -0.745 -0.198 0.081 mean(anti.table$Age)## [1] 2.226901anti.table$age.gc<-anti.table$Age-mean(anti.table$Age)
### Temperature
aggregate(temp~misuse,og,FUN=mean) fit.t <- glmer(misuse_category ~ WHO_participation*Temperature+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
fit.t.poly <- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
anova(fit.t,fit.t.poly) ## cubic term not required.
fit.t.3 <- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+I(Temperature**3)+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
anova(fit.t.3,fit.t.poly)summary(fit.t.poly)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## misuse_category ~ WHO_participation * Temperature + I(Temperature^2) +
## (1 | Hospital/Doctor.Identifier)
## Data: anti.table
## Control:
## glmerControl(optimizer = "bobyqa", boundary.tol = 3e-05, optCtrl = list(maxfun = 3e+05))
##
## AIC BIC logLik deviance df.resid
## 966.3 999.5 -476.1 952.3 848
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3900 -0.8222 0.4316 0.6060 1.7575
##
## Random effects:
## Groups Name Variance Std.Dev.
## Doctor.Identifier:Hospital (Intercept) 0.1161 0.3407
## Hospital (Intercept) 0.3036 0.5510
## Number of obs: 855, groups: Doctor.Identifier:Hospital, 134; Hospital, 36
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) 1.88238 0.30029 6.268
## WHO_participationWHO participation -1.64670 0.52249 -3.152
## Temperature -0.04521 0.31646 -0.143
## I(Temperature^2) -0.19073 0.08610 -2.215
## WHO_participationWHO participation:Temperature 0.41173 0.22798 1.806
## Pr(>|z|)
## (Intercept) 3.65e-10 ***
## WHO_participationWHO participation 0.00162 **
## Temperature 0.88641
## I(Temperature^2) 0.02674 *
## WHO_participationWHO participation:Temperature 0.07092 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WHO_pWHOp Tmprtr I(T^2)
## WHO_prtWHOp -0.293
## Temperature -0.763 0.007
## I(Tmprtr^2) 0.592 0.109 -0.946
## WHO_pWHOp:T 0.193 -0.789 -0.017 -0.130Level 1 Self Payment
- Patient self payment indicator was not significant as an independent fixed effect (p=0.59) nor as an interactive term with WHO participation status (p=0.86), whereas Temperature was a marginally significant interactive term (p=0.0719) in this joint model.
- The addition of Patient self payment increased the AIC compared to the temperature model (increase AIC is 3.71, Total AIC: 970), and was not included into the model.
##paymed
aggregate(Paymed~misuse,og,FUN=mean)fit.t.poly.paymed <- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+WHO_participation*Patient_Self_Payment+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
fit.t.paymed <- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+Patient_Self_Payment+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))Level 1 Self medication
- Prior self medication was significant (p=0.0011) as a fixed effect but not as an interactive term (p=0.98). Upon adding this to the temperature model, the AIC reduced, and the ANOVA suggested the models were different (p=0.002). the updated AIC 957.9, whereas the temerature model AIC was 966.3.
- Including Prior self medication into the temperature model modified the WHO effect estimate 12.6%, suggesting we should include this term into the further models.
aggregate(Selfmed~misuse,og,FUN=mean)fit.selfmed <- glmer(misuse_category ~ WHO_participation*Prior_self_medication+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
## includding with temp
fit.t.poly.selfmed <- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+WHO_participation*Prior_self_medication+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
anova(fit.t.poly,fit.t.poly.selfmed) (-1.64670 +1.85539)/(-1.64670)## [1] -0.1267323Wrong diagnosis
- Wrong diagnosis was a significant factor added separately to the model (p<2e-16) and we included this into the multivariate model.
- Including wrong diagnosis improved the previous model multivariate model and reduced AIC dramatically (AIC: 853.8). Interactions between WHO and wrong diagnosis were not significant (p=0.82).
- Interactions between wrong diagnosis and patient self payment were also not observed, indicating that the diagnosis is not modified by the patient payment method (p=0.23).
## WR diag
aggregate(Wrdiag~misuse,og,FUN=mean) fit.wrd <- glmer(misuse_category ~ WHO_participation*Wrong_diagnosis+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
## includding with temp
## AIC: 853.8
fit.t.self.wrd <- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+WHO_participation*Prior_self_medication+Wrong_diagnosis+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
## there was not an interaction between WHO and Wrong diagnosis.
fit.t.self.wrd2 <- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+WHO_participation*Prior_self_medication+Wrong_diagnosis*WHO_participation+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
anova(fit.t.self.wrd ,fit.t.self.wrd2) anova(fit.t.self.wrd ,fit.t.poly.selfmed) summary(fit.t.self.wrd2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## misuse_category ~ WHO_participation * Temperature + I(Temperature^2) +
## WHO_participation * Prior_self_medication + Wrong_diagnosis *
## WHO_participation + (1 | Hospital/Doctor.Identifier)
## Data: anti.table
## Control:
## glmerControl(optimizer = "bobyqa", boundary.tol = 3e-05, optCtrl = list(maxfun = 3e+05))
##
## AIC BIC logLik deviance df.resid
## 855.7 908.0 -416.9 833.7 844
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0690 -0.6063 0.3080 0.5374 2.3953
##
## Random effects:
## Groups Name Variance Std.Dev.
## Doctor.Identifier:Hospital (Intercept) 0.08404 0.2899
## Hospital (Intercept) 0.28749 0.5362
## Number of obs: 855, groups: Doctor.Identifier:Hospital, 134; Hospital, 36
##
## Fixed effects:
## Estimate
## (Intercept) 0.83440
## WHO_participationWHO participation -1.99439
## Temperature -0.01622
## I(Temperature^2) -0.21027
## Prior_self_medicationPatient self medicating -0.61409
## Wrong_diagnosisWrong diagnosis 1.94145
## WHO_participationWHO participation:Temperature 0.42960
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating 0.19538
## WHO_participationWHO participation:Wrong_diagnosisWrong diagnosis -0.14323
## Std. Error
## (Intercept) 0.34301
## WHO_participationWHO participation 0.71667
## Temperature 0.34207
## I(Temperature^2) 0.09247
## Prior_self_medicationPatient self medicating 0.22424
## Wrong_diagnosisWrong diagnosis 0.21340
## WHO_participationWHO participation:Temperature 0.24182
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating 0.59346
## WHO_participationWHO participation:Wrong_diagnosisWrong diagnosis 0.63143
## z value
## (Intercept) 2.433
## WHO_participationWHO participation -2.783
## Temperature -0.047
## I(Temperature^2) -2.274
## Prior_self_medicationPatient self medicating -2.739
## Wrong_diagnosisWrong diagnosis 9.098
## WHO_participationWHO participation:Temperature 1.777
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating 0.329
## WHO_participationWHO participation:Wrong_diagnosisWrong diagnosis -0.227
## Pr(>|z|)
## (Intercept) 0.01499
## WHO_participationWHO participation 0.00539
## Temperature 0.96218
## I(Temperature^2) 0.02297
## Prior_self_medicationPatient self medicating 0.00617
## Wrong_diagnosisWrong diagnosis < 2e-16
## WHO_participationWHO participation:Temperature 0.07565
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating 0.74199
## WHO_participationWHO participation:Wrong_diagnosisWrong diagnosis 0.82056
##
## (Intercept) *
## WHO_participationWHO participation **
## Temperature
## I(Temperature^2) *
## Prior_self_medicationPatient self medicating **
## Wrong_diagnosisWrong diagnosis ***
## WHO_participationWHO participation:Temperature .
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating
## WHO_participationWHO participation:Wrong_diagnosisWrong diagnosis
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WHO_pWHOp Tmprtr I(T^2) P__Psm Wrn_Wd WHO_WHOp: WHOpsm
## WHO_prtWHOp -0.281
## Temperature -0.730 0.028
## I(Tmprtr^2) 0.580 0.061 -0.945
## Prr_slf_Psm -0.217 0.123 -0.021 0.049
## Wrng_dgnsWd -0.269 0.106 0.003 -0.069 -0.020
## WHO_pWHOp:T 0.178 -0.531 -0.035 -0.115 -0.044 0.080
## WHO_WHOp:sm 0.075 -0.085 0.037 -0.063 -0.371 0.034 -0.052
## WHO_WHOp:Wd 0.110 -0.630 -0.024 0.039 0.017 -0.329 -0.113 -0.039 fit.t.self.wrd3 <- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+WHO_participation*Prior_self_medication+Wrong_diagnosis*Patient_Self_Payment+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
summary(fit.t.self.wrd3)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## misuse_category ~ WHO_participation * Temperature + I(Temperature^2) +
## WHO_participation * Prior_self_medication + Wrong_diagnosis *
## Patient_Self_Payment + (1 | Hospital/Doctor.Identifier)
## Data: anti.table
## Control:
## glmerControl(optimizer = "bobyqa", boundary.tol = 3e-05, optCtrl = list(maxfun = 3e+05))
##
## AIC BIC logLik deviance df.resid
## 855.3 912.3 -415.6 831.3 843
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0642 -0.5856 0.3118 0.5237 2.3900
##
## Random effects:
## Groups Name Variance Std.Dev.
## Doctor.Identifier:Hospital (Intercept) 0.07332 0.2708
## Hospital (Intercept) 0.32129 0.5668
## Number of obs: 855, groups: Doctor.Identifier:Hospital, 134; Hospital, 36
##
## Fixed effects:
## Estimate
## (Intercept) 0.90086
## WHO_participationWHO participation -2.08772
## Temperature -0.03790
## I(Temperature^2) -0.20836
## Prior_self_medicationPatient self medicating -0.64815
## Wrong_diagnosisWrong diagnosis 1.82462
## Patient_Self_PaymentPatient self paying -0.10271
## WHO_participationWHO participation:Temperature 0.44982
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating 0.03596
## Wrong_diagnosisWrong diagnosis:Patient_Self_PaymentPatient self paying 0.61575
## Std. Error
## (Intercept) 0.35350
## WHO_participationWHO participation 0.56330
## Temperature 0.34527
## I(Temperature^2) 0.09297
## Prior_self_medicationPatient self medicating 0.22733
## Wrong_diagnosisWrong diagnosis 0.21929
## Patient_Self_PaymentPatient self paying 0.45351
## WHO_participationWHO participation:Temperature 0.24223
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating 0.60815
## Wrong_diagnosisWrong diagnosis:Patient_Self_PaymentPatient self paying 0.51092
## z value
## (Intercept) 2.548
## WHO_participationWHO participation -3.706
## Temperature -0.110
## I(Temperature^2) -2.241
## Prior_self_medicationPatient self medicating -2.851
## Wrong_diagnosisWrong diagnosis 8.320
## Patient_Self_PaymentPatient self paying -0.226
## WHO_participationWHO participation:Temperature 1.857
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating 0.059
## Wrong_diagnosisWrong diagnosis:Patient_Self_PaymentPatient self paying 1.205
## Pr(>|z|)
## (Intercept) 0.01082
## WHO_participationWHO participation 0.00021
## Temperature 0.91260
## I(Temperature^2) 0.02502
## Prior_self_medicationPatient self medicating 0.00436
## Wrong_diagnosisWrong diagnosis < 2e-16
## Patient_Self_PaymentPatient self paying 0.82083
## WHO_participationWHO participation:Temperature 0.06331
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating 0.95285
## Wrong_diagnosisWrong diagnosis:Patient_Self_PaymentPatient self paying 0.22813
##
## (Intercept) *
## WHO_participationWHO participation ***
## Temperature
## I(Temperature^2) *
## Prior_self_medicationPatient self medicating **
## Wrong_diagnosisWrong diagnosis ***
## Patient_Self_PaymentPatient self paying
## WHO_participationWHO participation:Temperature .
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating
## Wrong_diagnosisWrong diagnosis:Patient_Self_PaymentPatient self paying
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WHO_pWHOp Tmprtr I(T^2) P__Psm Wrn_Wd P_S_sp WHO_WHOp:
## WHO_prtWHOp -0.280
## Temperature -0.715 0.020
## I(Tmprtr^2) 0.565 0.105 -0.945
## Prr_slf_Psm -0.211 0.164 -0.012 0.041
## Wrng_dgnsWd -0.304 -0.106 -0.003 -0.045 -0.006
## Ptnt_S_PPsp -0.196 0.054 -0.010 0.018 -0.053 0.323
## WHO_pWHOp:T 0.202 -0.773 -0.045 -0.105 -0.046 0.012 -0.038
## WHO_WHOp:sm 0.071 -0.137 0.040 -0.059 -0.340 0.043 -0.035 -0.068
## W_Wd:P_S_sp 0.203 -0.036 -0.020 -0.010 -0.033 -0.385 -0.741 0.074
## WHOpsm
## WHO_prtWHOp
## Temperature
## I(Tmprtr^2)
## Prr_slf_Psm
## Wrng_dgnsWd
## Ptnt_S_PPsp
## WHO_pWHOp:T
## WHO_WHOp:sm
## W_Wd:P_S_sp -0.089 ################Doctor education
- Doctor education level fitted into a model that included WHO and education, was not significant and did not differ from the variance component model (p=0.86).
- random slope on WHO could not converge.
- Doctor education levels were not significant when included into the multivariate model (AIC:919).
- including interaction terms between wrong diagnosis and doctor education levels struggled with convergence, but also had higher AIC approximately 855-860, whereas the previous MLM had 853.8.
- we did not detect interactions between WHO participation and doctor education levels (p>0.05) nor interactions between wrong diagnosis and education (p>0.05).
- After adding the doctor education level, the effect modification on the WHO estimate was approximately 8%, not enough to conclude confounding.
- Including doctor education pooled categories did not have significant associations compared against low-education, but the AIC was approximately equal to the model without education (AIC: 853.6). ANOVA p-value=0.0416, suggesting a marginal improvement.
##LEVEL 2
aggregate(DRed~misuse,og,FUN=mean) fit.dred <- glmer(misuse_category ~ WHO_participation*Doctor_education_Level+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
anova(fit3,fit.dred) ##pooling variables.
anti.table$Doctor_education_Level2<-as.character(anti.table$Doctor_education_Level)
anti.table[which(anti.table$Doctor_education_Level2%in%c("2","1")),"Doctor_education_Level2"]<-'Low education'
anti.table[which(anti.table$Doctor_education_Level2%in%c("4","3")),"Doctor_education_Level2"]<-'moderate education'
anti.table[which(anti.table$Doctor_education_Level2%in%c("5","6")),"Doctor_education_Level2"]<-'Highest education'
table(anti.table$Doctor_education_Level,anti.table$Doctor_education_Level2)##
## Highest education Low education moderate education
## 1 0 11 0
## 2 0 250 0
## 3 0 0 73
## 4 0 0 394
## 5 39 0 0
## 6 88 0 0 anti.table$Doctor_education_Level2<-factor(anti.table$Doctor_education_Level2)
anti.table$Doctor_education_Level2<-relevel(anti.table$Doctor_education_Level2,ref='Low education')
#education labels not significant.
fit.dred2 <- glmer(misuse_category ~ WHO_participation*Doctor_education_Level2+ Wrong_diagnosis*Doctor_education_Level2+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
### wrong diag interaction with DR ed
##AIC 853
fit.wd.ed <- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+WHO_participation*Prior_self_medication+Wrong_diagnosis+Doctor_education_Level2+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
anova(fit.t.self.wrd2,fit.wd.ed) ## AIC: 857.3
fit.wd.ed.inter <- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+WHO_participation*Prior_self_medication+Wrong_diagnosis*Doctor_education_Level2+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
# adding random slope WHO (would not converge).
##
## adding random slope of DR ED.
# fit.t.self.wrd.slope <- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+WHO_participation*Prior_self_medication+Wrong_diagnosis+(Doctor_education_Level2|Hospital/Doctor.Identifier),
# data = anti.table,
# family = binomial,
# control=glmerControl(optimizer ="Nelder_Mead",boundary.tol=3e-7,
# optCtrl = list(maxfun=3e6)))
aggregate(WHO~misuse,og,FUN=mean)Re-examining Age at the final model
- The anova comparison of the simple addition of age, compared to the variance component nested model (p=0.18) did not identify age as an important feature.
- However adding in age (grand mean centered) feature into the multivariate model was significant (p=0.02).
- Age did not interact with wrong diagnosis (p=0.19)
- Age was not a quadratic polynomial feature (p=0.63).
- The AIC of the model that includes temperature, quadratic temperature (p=0.02), self medication (0.004), wrong diagnosis (p<2e-16), age (p=0.03) and doctor education is 849.3.
anova(fit.age,fit3) ##interaction with wrong diag and age
fit.wd.age.inter <- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+WHO_participation*Prior_self_medication+Wrong_diagnosis*age.gc+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
summary(fit.wd.age.inter)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## misuse_category ~ WHO_participation * Temperature + I(Temperature^2) +
## WHO_participation * Prior_self_medication + Wrong_diagnosis *
## age.gc + (1 | Hospital/Doctor.Identifier)
## Data: anti.table
## Control:
## glmerControl(optimizer = "bobyqa", boundary.tol = 3e-05, optCtrl = list(maxfun = 3e+05))
##
## AIC BIC logLik deviance df.resid
## 851.1 908.1 -413.6 827.1 843
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0202 -0.5982 0.3046 0.5201 2.6834
##
## Random effects:
## Groups Name Variance Std.Dev.
## Doctor.Identifier:Hospital (Intercept) 0.07816 0.2796
## Hospital (Intercept) 0.27826 0.5275
## Number of obs: 855, groups: Doctor.Identifier:Hospital, 134; Hospital, 36
##
## Fixed effects:
## Estimate
## (Intercept) 0.846936
## WHO_participationWHO participation -2.051018
## Temperature -0.008698
## I(Temperature^2) -0.213126
## Prior_self_medicationPatient self medicating -0.660340
## Wrong_diagnosisWrong diagnosis 1.921005
## age.gc 0.274766
## WHO_participationWHO participation:Temperature 0.429004
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating 0.134868
## Wrong_diagnosisWrong diagnosis:age.gc -0.185307
## Std. Error
## (Intercept) 0.341188
## WHO_participationWHO participation 0.559024
## Temperature 0.343612
## I(Temperature^2) 0.093039
## Prior_self_medicationPatient self medicating 0.224769
## Wrong_diagnosisWrong diagnosis 0.202947
## age.gc 0.115960
## WHO_participationWHO participation:Temperature 0.241899
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating 0.599894
## Wrong_diagnosisWrong diagnosis:age.gc 0.143418
## z value
## (Intercept) 2.482
## WHO_participationWHO participation -3.669
## Temperature -0.025
## I(Temperature^2) -2.291
## Prior_self_medicationPatient self medicating -2.938
## Wrong_diagnosisWrong diagnosis 9.466
## age.gc 2.369
## WHO_participationWHO participation:Temperature 1.773
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating 0.225
## Wrong_diagnosisWrong diagnosis:age.gc -1.292
## Pr(>|z|)
## (Intercept) 0.013053
## WHO_participationWHO participation 0.000244
## Temperature 0.979805
## I(Temperature^2) 0.021980
## Prior_self_medicationPatient self medicating 0.003305
## Wrong_diagnosisWrong diagnosis < 2e-16
## age.gc 0.017813
## WHO_participationWHO participation:Temperature 0.076148
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating 0.822120
## Wrong_diagnosisWrong diagnosis:age.gc 0.196331
##
## (Intercept) *
## WHO_participationWHO participation ***
## Temperature
## I(Temperature^2) *
## Prior_self_medicationPatient self medicating **
## Wrong_diagnosisWrong diagnosis ***
## age.gc *
## WHO_participationWHO participation:Temperature .
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating
## Wrong_diagnosisWrong diagnosis:age.gc
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WHO_pWHOp Tmprtr I(T^2) P__Psm Wrn_Wd age.gc WHO_WHOp:
## WHO_prtWHOp -0.267
## Temperature -0.730 0.013
## I(Tmprtr^2) 0.579 0.111 -0.946
## Prr_slf_Psm -0.223 0.167 -0.022 0.052
## Wrng_dgnsWd -0.261 -0.132 0.000 -0.062 -0.015
## age.gc 0.018 -0.018 0.007 -0.024 -0.077 0.029
## WHO_pWHOp:T 0.191 -0.780 -0.036 -0.110 -0.047 0.041 0.017
## WHO_WHOp:sm 0.076 -0.139 0.039 -0.063 -0.357 0.016 0.025 -0.070
## Wrng_dgWd:. -0.010 0.051 -0.007 0.013 0.020 -0.007 -0.791 0.000
## WHOpsm
## WHO_prtWHOp
## Temperature
## I(Tmprtr^2)
## Prr_slf_Psm
## Wrng_dgnsWd
## age.gc
## WHO_pWHOp:T
## WHO_WHOp:sm
## Wrng_dgWd:. -0.101##age as a quadratic term is not significant.
fit.wd.age.add <- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+age.gc+I(age.gc**2)+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
summary(fit.wd.age.add)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## misuse_category ~ WHO_participation * Temperature + I(Temperature^2) +
## Prior_self_medication + Wrong_diagnosis + age.gc + I(age.gc^2) +
## (1 | Hospital/Doctor.Identifier)
## Data: anti.table
## Control:
## glmerControl(optimizer = "bobyqa", boundary.tol = 3e-05, optCtrl = list(maxfun = 3e+05))
##
## AIC BIC logLik deviance df.resid
## 850.6 902.9 -414.3 828.6 844
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1049 -0.6038 0.3100 0.5361 2.5067
##
## Random effects:
## Groups Name Variance Std.Dev.
## Doctor.Identifier:Hospital (Intercept) 0.0802 0.2832
## Hospital (Intercept) 0.2558 0.5057
## Number of obs: 855, groups: Doctor.Identifier:Hospital, 134; Hospital, 36
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) 0.89867 0.35739 2.515
## WHO_participationWHO participation -2.01405 0.54979 -3.663
## Temperature -0.02033 0.34319 -0.059
## I(Temperature^2) -0.21012 0.09290 -2.262
## Prior_self_medicationPatient self medicating -0.65650 0.21041 -3.120
## Wrong_diagnosisWrong diagnosis 1.92481 0.20205 9.526
## age.gc 0.14974 0.07198 2.080
## I(age.gc^2) -0.02680 0.05568 -0.481
## WHO_participationWHO participation:Temperature 0.42785 0.24186 1.769
## Pr(>|z|)
## (Intercept) 0.011918 *
## WHO_participationWHO participation 0.000249 ***
## Temperature 0.952751
## I(Temperature^2) 0.023702 *
## Prior_self_medicationPatient self medicating 0.001808 **
## Wrong_diagnosisWrong diagnosis < 2e-16 ***
## age.gc 0.037491 *
## I(age.gc^2) 0.630221
## WHO_participationWHO participation:Temperature 0.076894 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WHO_pWHOp Tmprtr I(T^2) P__Psm Wrn_Wd age.gc I(.^2)
## WHO_prtWHOp -0.244
## Temperature -0.713 0.020
## I(Tmprtr^2) 0.568 0.102 -0.946
## Prr_slf_Psm -0.226 0.127 -0.003 0.026
## Wrng_dgnsWd -0.244 -0.132 0.000 -0.062 -0.014
## age.gc -0.052 0.024 0.017 -0.035 -0.117 0.051
## I(age.gc^2) -0.327 0.000 0.042 -0.031 0.078 -0.001 0.238
## WHO_pWHOp:T 0.175 -0.804 -0.033 -0.114 -0.071 0.042 0.029 0.038## including DR education improves the model.
fit.wd.age.add.ed <- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+age.gc+Doctor_education_Level2+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
summary(fit.wd.age.add.ed)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## misuse_category ~ WHO_participation * Temperature + I(Temperature^2) +
## Prior_self_medication + Wrong_diagnosis + age.gc + Doctor_education_Level2 +
## (1 | Hospital/Doctor.Identifier)
## Data: anti.table
## Control:
## glmerControl(optimizer = "bobyqa", boundary.tol = 3e-05, optCtrl = list(maxfun = 3e+05))
##
## AIC BIC logLik deviance df.resid
## 849.3 906.3 -412.6 825.3 843
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1975 -0.5864 0.3049 0.5386 2.5563
##
## Random effects:
## Groups Name Variance Std.Dev.
## Doctor.Identifier:Hospital (Intercept) 0.04089 0.2022
## Hospital (Intercept) 0.27254 0.5221
## Number of obs: 855, groups: Doctor.Identifier:Hospital, 134; Hospital, 36
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) 0.85206 0.38130 2.235
## WHO_participationWHO participation -1.87258 0.55209 -3.392
## Temperature -0.01565 0.34076 -0.046
## I(Temperature^2) -0.21233 0.09228 -2.301
## Prior_self_medicationPatient self medicating -0.61396 0.20889 -2.939
## Wrong_diagnosisWrong diagnosis 1.95478 0.20274 9.642
## age.gc 0.14828 0.07008 2.116
## Doctor_education_Level2Highest education -0.54454 0.35890 -1.517
## Doctor_education_Level2moderate education 0.03361 0.26461 0.127
## WHO_participationWHO participation:Temperature 0.39459 0.24132 1.635
## Pr(>|z|)
## (Intercept) 0.025441 *
## WHO_participationWHO participation 0.000694 ***
## Temperature 0.963370
## I(Temperature^2) 0.021397 *
## Prior_self_medicationPatient self medicating 0.003291 **
## Wrong_diagnosisWrong diagnosis < 2e-16 ***
## age.gc 0.034359 *
## Doctor_education_Level2Highest education 0.129196
## Doctor_education_Level2moderate education 0.898931
## WHO_participationWHO participation:Temperature 0.102014
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WHO_pWHOp Tmprtr I(T^2) P__Psm Wrn_Wd age.gc D__L2He Dc__L2e
## WHO_prtWHOp -0.241
## Temperature -0.658 0.018
## I(Tmprtr^2) 0.516 0.102 -0.945
## Prr_slf_Psm -0.178 0.131 -0.008 0.030
## Wrng_dgnsWd -0.182 -0.112 -0.003 -0.063 -0.001
## age.gc 0.024 0.014 0.008 -0.030 -0.143 0.046
## Dctr_d_L2He -0.306 -0.076 0.016 0.016 -0.073 -0.161 0.047
## Dctr_dc_L2e -0.465 0.035 0.014 0.004 -0.016 -0.076 -0.003 0.550
## WHO_pWHOp:T 0.158 -0.802 -0.031 -0.117 -0.081 0.028 0.027 0.077 0.028Self Payment interactions
- we hypothesized that self-payment may be interaction term with wrong diagnosis or doctor level of education, however it was not.
- The AIC of the multi-level payment feature was 918.
- We did not find interactions with doctor level features and payment method.
## Payment interactions with education level
fit.pay.inters <- glmer(misuse_category ~ WHO_participation+Wrong_diagnosis*Patient_Self_Payment+Patient_Self_Payment*Doctor_education_Level2+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
summary(fit.pay.inters)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## misuse_category ~ WHO_participation + Wrong_diagnosis * Patient_Self_Payment +
## Patient_Self_Payment * Doctor_education_Level2 + (1 | Hospital/Doctor.Identifier)
## Data: anti.table
## Control:
## glmerControl(optimizer = "bobyqa", boundary.tol = 3e-05, optCtrl = list(maxfun = 3e+05))
##
## AIC BIC logLik deviance df.resid
## 918.4 970.7 -448.2 896.4 844
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1814 -0.7554 0.3895 0.4877 3.0766
##
## Random effects:
## Groups Name Variance Std.Dev.
## Doctor.Identifier:Hospital (Intercept) 0.06904 0.2627
## Hospital (Intercept) 0.22243 0.4716
## Number of obs: 855, groups: Doctor.Identifier:Hospital, 134; Hospital, 36
##
## Fixed effects:
## Estimate
## (Intercept) -0.2411
## WHO_participationWHO participation -1.2378
## Wrong_diagnosisWrong diagnosis 1.8644
## Patient_Self_PaymentPatient self paying 1.6263
## Doctor_education_Level2Highest education -0.3594
## Doctor_education_Level2moderate education 0.1616
## Wrong_diagnosisWrong diagnosis:Patient_Self_PaymentPatient self paying 0.2838
## Patient_Self_PaymentPatient self paying:Doctor_education_Level2Highest education -1.6979
## Patient_Self_PaymentPatient self paying:Doctor_education_Level2moderate education -1.6271
## Std. Error
## (Intercept) 0.2405
## WHO_participationWHO participation 0.3105
## Wrong_diagnosisWrong diagnosis 0.2070
## Patient_Self_PaymentPatient self paying 1.4258
## Doctor_education_Level2Highest education 0.3786
## Doctor_education_Level2moderate education 0.2455
## Wrong_diagnosisWrong diagnosis:Patient_Self_PaymentPatient self paying 0.4924
## Patient_Self_PaymentPatient self paying:Doctor_education_Level2Highest education 1.4994
## Patient_Self_PaymentPatient self paying:Doctor_education_Level2moderate education 1.4556
## z value
## (Intercept) -1.003
## WHO_participationWHO participation -3.986
## Wrong_diagnosisWrong diagnosis 9.007
## Patient_Self_PaymentPatient self paying 1.141
## Doctor_education_Level2Highest education -0.949
## Doctor_education_Level2moderate education 0.658
## Wrong_diagnosisWrong diagnosis:Patient_Self_PaymentPatient self paying 0.576
## Patient_Self_PaymentPatient self paying:Doctor_education_Level2Highest education -1.132
## Patient_Self_PaymentPatient self paying:Doctor_education_Level2moderate education -1.118
## Pr(>|z|)
## (Intercept) 0.316
## WHO_participationWHO participation 6.72e-05
## Wrong_diagnosisWrong diagnosis < 2e-16
## Patient_Self_PaymentPatient self paying 0.254
## Doctor_education_Level2Highest education 0.342
## Doctor_education_Level2moderate education 0.510
## Wrong_diagnosisWrong diagnosis:Patient_Self_PaymentPatient self paying 0.564
## Patient_Self_PaymentPatient self paying:Doctor_education_Level2Highest education 0.257
## Patient_Self_PaymentPatient self paying:Doctor_education_Level2moderate education 0.264
##
## (Intercept)
## WHO_participationWHO participation ***
## Wrong_diagnosisWrong diagnosis ***
## Patient_Self_PaymentPatient self paying
## Doctor_education_Level2Highest education
## Doctor_education_Level2moderate education
## Wrong_diagnosisWrong diagnosis:Patient_Self_PaymentPatient self paying
## Patient_Self_PaymentPatient self paying:Doctor_education_Level2Highest education
## Patient_Self_PaymentPatient self paying:Doctor_education_Level2moderate education
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WHO_Wp Wrn_Wd P_S_sp D__L2He Dc__L2e W_Wdsp
## WHO_prtWHOp -0.238
## Wrng_dgnsWd -0.444 -0.154
## Ptnt_S_PPsp -0.144 -0.022 0.101
## Dctr_d_L2He -0.361 -0.028 -0.134 0.073
## Dctr_dc_L2e -0.625 0.072 -0.060 0.114 0.480
## W_Wd:P_S_sp 0.189 0.032 -0.397 -0.112 0.045 0.024
## P_S_PPsp:D__L2He 0.074 0.031 0.018 -0.918 -0.210 -0.112 -0.151
## Pt_S_PPsp:D__L2e 0.082 0.028 -0.003 -0.949 -0.067 -0.151 -0.094
## P_S_PPsp:D__L2He
## WHO_prtWHOp
## Wrng_dgnsWd
## Ptnt_S_PPsp
## Dctr_d_L2He
## Dctr_dc_L2e
## W_Wd:P_S_sp
## P_S_PPsp:D__L2He
## Pt_S_PPsp:D__L2e 0.926## Payment interactions with education level
fit.pay.inters2 <- glmer(misuse_category ~ WHO_participation+Wrong_diagnosis*Patient_Self_Payment+Patient_Self_Payment*Wrong_diagnosis+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
summary(fit.pay.inters2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## misuse_category ~ WHO_participation + Wrong_diagnosis * Patient_Self_Payment +
## Patient_Self_Payment * Wrong_diagnosis + (1 | Hospital/Doctor.Identifier)
## Data: anti.table
## Control:
## glmerControl(optimizer = "bobyqa", boundary.tol = 3e-05, optCtrl = list(maxfun = 3e+05))
##
## AIC BIC logLik deviance df.resid
## 915.5 948.8 -450.8 901.5 848
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0783 -0.7896 0.4001 0.4776 2.6361
##
## Random effects:
## Groups Name Variance Std.Dev.
## Doctor.Identifier:Hospital (Intercept) 0.1007 0.3173
## Hospital (Intercept) 0.2186 0.4675
## Number of obs: 855, groups: Doctor.Identifier:Hospital, 134; Hospital, 36
##
## Fixed effects:
## Estimate
## (Intercept) -0.17308
## WHO_participationWHO participation -1.30000
## Wrong_diagnosisWrong diagnosis 1.85870
## Patient_Self_PaymentPatient self paying 0.03977
## Wrong_diagnosisWrong diagnosis:Patient_Self_PaymentPatient self paying 0.18067
## Std. Error
## (Intercept) 0.18798
## WHO_participationWHO participation 0.31022
## Wrong_diagnosisWrong diagnosis 0.20675
## Patient_Self_PaymentPatient self paying 0.42592
## Wrong_diagnosisWrong diagnosis:Patient_Self_PaymentPatient self paying 0.48062
## z value
## (Intercept) -0.921
## WHO_participationWHO participation -4.191
## Wrong_diagnosisWrong diagnosis 8.990
## Patient_Self_PaymentPatient self paying 0.093
## Wrong_diagnosisWrong diagnosis:Patient_Self_PaymentPatient self paying 0.376
## Pr(>|z|)
## (Intercept) 0.357
## WHO_participationWHO participation 2.78e-05
## Wrong_diagnosisWrong diagnosis < 2e-16
## Patient_Self_PaymentPatient self paying 0.926
## Wrong_diagnosisWrong diagnosis:Patient_Self_PaymentPatient self paying 0.707
##
## (Intercept)
## WHO_participationWHO participation ***
## Wrong_diagnosisWrong diagnosis ***
## Patient_Self_PaymentPatient self paying
## Wrong_diagnosisWrong diagnosis:Patient_Self_PaymentPatient self paying
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WHO_Wp Wrn_Wd P_S_sp
## WHO_prtWHOp -0.250
## Wrng_dgnsWd -0.640 -0.161
## Ptnt_S_PPsp -0.361 0.028 0.343
## W_Wd:P_S_sp 0.277 0.027 -0.403 -0.758Random slopes
- The slopes fitted with WHO participation would not converge.
- adding doctor level education would not improve the model, and would increase the AIC if we added education. Education as a random slope would not converge.
- including wrong diagnosis as a random slope, improves the model (AIC: 844.9), and the anova is significantly different from the random intercept model only (p=0.007). There is a correlation between hospitals and patient diagnosis.
- The interaction between patient self medication and WHO participant was not significant (p=0.74), and keeping the self medication as an additive effect improved the model and reduced the AIC : 843.0 ( and the anova comparison had p=0.7466).
- The main effect of temperature is not significant (p=0.9), however the interaction between WHO and temperature is significant (p<0.05),if we only keep the interaction term and drop the main effect, temperature/WHO interaction is not significant (p=0.3). Hence we keep the main effect of temperature.
##
dat<-NULL
for(h in unique(anti.table$Hospital)){
x<-anti.table[which(anti.table$Hospital==h),]
if(nrow(x)>=4){
cc<-as.numeric(coef(glm(misuse_category~Wrong_diagnosis,family='binomial',data=x)))
cc<-c(cc,h)
dat<-rbind(dat,cc)
}
}
dat[,1]<-as.numeric(dat[,1])
dat[,2]<-as.numeric(dat[,2])
dat<-as.data.frame(dat)
colnames(dat)<-c("Int","slopes","Hospital")
ggplot(dat,aes(x=Int,y=slopes))+geom_point()+geom_smooth()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
slope<- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+WHO_participation*Prior_self_medication+Wrong_diagnosis+Age+(Wrong_diagnosis|Hospital) +(1|Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
summary(slope)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## misuse_category ~ WHO_participation * Temperature + I(Temperature^2) +
## WHO_participation * Prior_self_medication + Wrong_diagnosis +
## Age + (Wrong_diagnosis | Hospital) + (1 | Doctor.Identifier)
## Data: anti.table
## Control:
## glmerControl(optimizer = "bobyqa", boundary.tol = 3e-05, optCtrl = list(maxfun = 3e+05))
##
## AIC BIC logLik deviance df.resid
## 844.9 906.7 -409.5 818.9 842
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.4155 -0.5618 0.2941 0.5256 3.2537
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Doctor.Identifier (Intercept) 0.0372 0.1929
## Hospital (Intercept) 0.4966 0.7047
## Wrong_diagnosisWrong diagnosis 0.9643 0.9820 -0.69
## Number of obs: 855, groups: Doctor.Identifier, 134; Hospital, 36
##
## Fixed effects:
## Estimate
## (Intercept) 0.43725
## WHO_participationWHO participation -2.02635
## Temperature 0.02034
## I(Temperature^2) -0.22518
## Prior_self_medicationPatient self medicating -0.76248
## Wrong_diagnosisWrong diagnosis 2.04509
## Age 0.15935
## WHO_participationWHO participation:Temperature 0.41996
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating 0.19407
## Std. Error
## (Intercept) 0.38846
## WHO_participationWHO participation 0.57508
## Temperature 0.35007
## I(Temperature^2) 0.09508
## Prior_self_medicationPatient self medicating 0.23276
## Wrong_diagnosisWrong diagnosis 0.28076
## Age 0.07183
## WHO_participationWHO participation:Temperature 0.24758
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating 0.60006
## z value
## (Intercept) 1.126
## WHO_participationWHO participation -3.524
## Temperature 0.058
## I(Temperature^2) -2.368
## Prior_self_medicationPatient self medicating -3.276
## Wrong_diagnosisWrong diagnosis 7.284
## Age 2.218
## WHO_participationWHO participation:Temperature 1.696
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating 0.323
## Pr(>|z|)
## (Intercept) 0.260325
## WHO_participationWHO participation 0.000426
## Temperature 0.953668
## I(Temperature^2) 0.017865
## Prior_self_medicationPatient self medicating 0.001053
## Wrong_diagnosisWrong diagnosis 3.24e-13
## Age 0.026540
## WHO_participationWHO participation:Temperature 0.089842
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating 0.746372
##
## (Intercept)
## WHO_participationWHO participation ***
## Temperature
## I(Temperature^2) *
## Prior_self_medicationPatient self medicating **
## Wrong_diagnosisWrong diagnosis ***
## Age *
## WHO_participationWHO participation:Temperature .
## WHO_participationWHO participation:Prior_self_medicationPatient self medicating
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WHO_pWHOp Tmprtr I(T^2) P__Psm Wrn_Wd Age WHO_WHOp:
## WHO_prtWHOp -0.234
## Temperature -0.649 0.015
## I(Tmprtr^2) 0.526 0.109 -0.945
## Prr_slf_Psm -0.157 0.166 -0.028 0.060
## Wrng_dgnsWd -0.348 -0.130 0.001 -0.052 -0.036
## Age -0.391 0.033 0.004 -0.027 -0.092 0.036
## WHO_pWHOp:T 0.158 -0.774 -0.035 -0.114 -0.043 0.024 0.033
## WHO_WHOp:sm 0.100 -0.139 0.045 -0.067 -0.373 0.025 -0.090 -0.065slope2<- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+Age+(Wrong_diagnosis|Hospital) +(1|Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
summary(slope2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## misuse_category ~ WHO_participation * Temperature + I(Temperature^2) +
## Prior_self_medication + Wrong_diagnosis + Age + (Wrong_diagnosis |
## Hospital) + (1 | Doctor.Identifier)
## Data: anti.table
## Control:
## glmerControl(optimizer = "bobyqa", boundary.tol = 3e-05, optCtrl = list(maxfun = 3e+05))
##
## AIC BIC logLik deviance df.resid
## 843.0 900.1 -409.5 819.0 843
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3812 -0.5649 0.2957 0.5256 3.2025
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Doctor.Identifier (Intercept) 0.03439 0.1855
## Hospital (Intercept) 0.49367 0.7026
## Wrong_diagnosisWrong diagnosis 0.95372 0.9766 -0.69
## Number of obs: 855, groups: Doctor.Identifier, 134; Hospital, 36
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) 0.42459 0.38588 1.100
## WHO_participationWHO participation -2.00052 0.56868 -3.518
## Temperature 0.01529 0.34948 0.044
## I(Temperature^2) -0.22315 0.09481 -2.354
## Prior_self_medicationPatient self medicating -0.73452 0.21613 -3.399
## Wrong_diagnosisWrong diagnosis 2.04317 0.27989 7.300
## Age 0.16148 0.07149 2.259
## WHO_participationWHO participation:Temperature 0.42495 0.24730 1.718
## Pr(>|z|)
## (Intercept) 0.271191
## WHO_participationWHO participation 0.000435 ***
## Temperature 0.965108
## I(Temperature^2) 0.018583 *
## Prior_self_medicationPatient self medicating 0.000677 ***
## Wrong_diagnosisWrong diagnosis 2.88e-13 ***
## Age 0.023907 *
## WHO_participationWHO participation:Temperature 0.085734 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WHO_pWHOp Tmprtr I(T^2) P__Psm Wrn_Wd Age
## WHO_prtWHOp -0.221
## Temperature -0.658 0.022
## I(Tmprtr^2) 0.537 0.101 -0.945
## Prr_slf_Psm -0.128 0.122 -0.014 0.039
## Wrng_dgnsWd -0.354 -0.127 0.000 -0.050 -0.027
## Age -0.386 0.019 0.008 -0.033 -0.135 0.039
## WHO_pWHOp:T 0.163 -0.793 -0.032 -0.118 -0.069 0.026 0.029## the temperature interaction requires the main effect of temperature.
slope3<- glmer(misuse_category ~WHO_participation+WHO_participation:Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+Age+(Wrong_diagnosis|Hospital) +(1|Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))Final model without random slopes.
- The model is defined as i,j,k represents pediatric patient, doctors, and hospital. and \(Y_{ik}\) is an indicator for 1 for misuse of at least one antibiotic, or 0 correct use.
- WRITE THE COMPLETE FORMULA HELP!!!!
- \(logit(Pr(Y_{ik}==1))= (\beta_{0}+U_{0k}+U_{0jk})+\beta_{1}WHO_k+ \beta_2Temp_{ijk}+\beta_3Temp^{2}_{ijk}+\beta_4SelfMed_{ijk}+\beta_5*WrongDiag_{ijk}+\beta_6(Age_{ijk}-2.23)+\beta_7WHO_{k}Temp_{ijk}\)
## AIC: 853.8
fit.t.self.wrd <- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+age.gc+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
##interaction with wrong diag and age
final<- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+age.gc+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
summary(final)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## misuse_category ~ WHO_participation * Temperature + I(Temperature^2) +
## Prior_self_medication + Wrong_diagnosis + age.gc + (1 | Hospital/Doctor.Identifier)
## Data: anti.table
## Control:
## glmerControl(optimizer = "bobyqa", boundary.tol = 3e-05, optCtrl = list(maxfun = 3e+05))
##
## AIC BIC logLik deviance df.resid
## 848.8 896.3 -414.4 828.8 845
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0716 -0.6036 0.3090 0.5339 2.5049
##
## Random effects:
## Groups Name Variance Std.Dev.
## Doctor.Identifier:Hospital (Intercept) 0.08423 0.2902
## Hospital (Intercept) 0.26070 0.5106
## Number of obs: 855, groups: Doctor.Identifier:Hospital, 134; Hospital, 36
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) 0.84289 0.33836 2.491
## WHO_participationWHO participation -2.01517 0.55049 -3.661
## Temperature -0.01340 0.34307 -0.039
## I(Temperature^2) -0.21163 0.09288 -2.279
## Prior_self_medicationPatient self medicating -0.64884 0.21014 -3.088
## Wrong_diagnosisWrong diagnosis 1.92578 0.20234 9.518
## age.gc 0.15808 0.07027 2.250
## WHO_participationWHO participation:Temperature 0.43249 0.24156 1.790
## Pr(>|z|)
## (Intercept) 0.012735 *
## WHO_participationWHO participation 0.000252 ***
## Temperature 0.968844
## I(Temperature^2) 0.022689 *
## Prior_self_medicationPatient self medicating 0.002018 **
## Wrong_diagnosisWrong diagnosis < 2e-16 ***
## age.gc 0.024472 *
## WHO_participationWHO participation:Temperature 0.073387 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WHO_pWHOp Tmprtr I(T^2) P__Psm Wrn_Wd age.gc
## WHO_prtWHOp -0.259
## Temperature -0.739 0.020
## I(Tmprtr^2) 0.590 0.102 -0.946
## Prr_slf_Psm -0.212 0.127 -0.006 0.029
## Wrng_dgnsWd -0.258 -0.132 -0.001 -0.062 -0.013
## age.gc 0.027 0.024 0.008 -0.030 -0.140 0.052
## WHO_pWHOp:T 0.198 -0.803 -0.035 -0.114 -0.074 0.043 0.020anova(final,slope)Evaluating the model residuals
Overall the residuals of WHO_participation, self medication and wrong diagnosis are moderately balanced around 0. Wrong diagnosis is skewed and shows signs of skewness. The age term is very balanced around 0, and the Temperature residuals shows symptoms of non-linearity, which is
anti.table$predicted<-predict(final,type='response')
anti.table$residuals<-residuals(final,type='response')
ggplot(anti.table,aes(x=WHO_participation,y=residuals))+geom_boxplot()
ggplot(anti.table,aes(x=Prior_self_medication,y=residuals))+geom_boxplot()
ggplot(anti.table,aes(x=Wrong_diagnosis,y=residuals))+geom_boxplot()
### age.gc is balanced around 0
ggplot(anti.table,aes(x=age.gc,y=residuals))+geom_point()
## in the Temperature main effect, we do observe non-linearity in the residuals which is why we included the quadratic term.
ggplot(anti.table,aes(x=Temperature,y=residuals))+geom_point()
Hospital
- The performance of the hospital each hospital crosses 0 in their intervals, so the misuse
library(lattice)
myrandef<-ranef(final,condVar=TRUE)
dotplot(myrandef)## $`Doctor.Identifier:Hospital`
##
## $Hospital
final model with random slope
- Including the random slope of Wrong diagnosis had AIC of 843, without random slope was 848, which is quite close to not including the random slope.
final.slope<- glmer(misuse_category ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+age.gc+(Wrong_diagnosis|Hospital) +(1|Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
summary(final.slope)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## misuse_category ~ WHO_participation * Temperature + I(Temperature^2) +
## Prior_self_medication + Wrong_diagnosis + age.gc + (Wrong_diagnosis |
## Hospital) + (1 | Doctor.Identifier)
## Data: anti.table
## Control:
## glmerControl(optimizer = "bobyqa", boundary.tol = 3e-05, optCtrl = list(maxfun = 3e+05))
##
## AIC BIC logLik deviance df.resid
## 843.0 900.1 -409.5 819.0 843
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3812 -0.5649 0.2957 0.5256 3.2025
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Doctor.Identifier (Intercept) 0.0344 0.1855
## Hospital (Intercept) 0.4937 0.7026
## Wrong_diagnosisWrong diagnosis 0.9537 0.9766 -0.69
## Number of obs: 855, groups: Doctor.Identifier, 134; Hospital, 36
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) 0.78419 0.35619 2.202
## WHO_participationWHO participation -2.00052 0.56868 -3.518
## Temperature 0.01529 0.34948 0.044
## I(Temperature^2) -0.22315 0.09481 -2.354
## Prior_self_medicationPatient self medicating -0.73452 0.21613 -3.399
## Wrong_diagnosisWrong diagnosis 2.04318 0.27989 7.300
## age.gc 0.16148 0.07149 2.259
## WHO_participationWHO participation:Temperature 0.42495 0.24730 1.718
## Pr(>|z|)
## (Intercept) 0.027691 *
## WHO_participationWHO participation 0.000435 ***
## Temperature 0.965109
## I(Temperature^2) 0.018583 *
## Prior_self_medicationPatient self medicating 0.000677 ***
## Wrong_diagnosisWrong diagnosis 2.88e-13 ***
## age.gc 0.023907 *
## WHO_participationWHO participation:Temperature 0.085732 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WHO_pWHOp Tmprtr I(T^2) P__Psm Wrn_Wd age.gc
## WHO_prtWHOp -0.231
## Temperature -0.709 0.022
## I(Tmprtr^2) 0.567 0.101 -0.945
## Prr_slf_Psm -0.200 0.122 -0.014 0.039
## Wrng_dgnsWd -0.366 -0.127 0.000 -0.050 -0.027
## age.gc 0.029 0.019 0.008 -0.033 -0.135 0.039
## WHO_pWHOp:T 0.190 -0.793 -0.032 -0.118 -0.069 0.026 0.029#Secondary/Sensitivity analyses (Temperature) * Temperature (p=0.002), and Temperature\(^2\) (p=0.01) are significant with the second categorization.
anti.table$misuse<-anti$misuse
anti.table$misuse_category2<-ifelse(anti.table$Antibiotic_misuse=='misuse of several antibiotics',1,0)
table(anti.table$misuse,anti.table$misuse_category2)##
## 0 1
## correct use 261 0
## misuse of 1 antibiotic 380 0
## misuse of several antibiotics 0 214 fit.t <- glmer(misuse_category2 ~ WHO_participation*Temperature+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
fit.t.poly <- glmer(misuse_category2 ~ WHO_participation*Temperature+I(Temperature**2)+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
anova(fit.t,fit.t.poly)Secondary analysis Payment
- Patient self payment is not significant (p=0.14).
fit.t.paymed <- glmer(misuse_category2 ~ WHO_participation*Temperature+I(Temperature**2)+Patient_Self_Payment+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))Secondary self medication
- Self medication is significant (p=0.02)
## includding with temp
fit.t.poly.selfmed <- glmer(misuse_category2 ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
summary(fit.t.poly.selfmed)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## misuse_category2 ~ WHO_participation * Temperature + I(Temperature^2) +
## Prior_self_medication + (1 | Hospital/Doctor.Identifier)
## Data: anti.table
## Control:
## glmerControl(optimizer = "bobyqa", boundary.tol = 3e-05, optCtrl = list(maxfun = 3e+05))
##
## AIC BIC logLik deviance df.resid
## 842.8 880.8 -413.4 826.8 847
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6660 -0.5120 -0.2794 0.2268 3.0929
##
## Random effects:
## Groups Name Variance Std.Dev.
## Doctor.Identifier:Hospital (Intercept) 0.842 0.9176
## Hospital (Intercept) 1.011 1.0055
## Number of obs: 855, groups: Doctor.Identifier:Hospital, 134; Hospital, 36
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) -2.1693 0.4062 -5.341
## WHO_participationWHO participation -3.3746 1.2223 -2.761
## Temperature 1.1935 0.3915 3.049
## I(Temperature^2) -0.2669 0.1076 -2.481
## Prior_self_medicationPatient self medicating -0.6027 0.2522 -2.390
## WHO_participationWHO participation:Temperature 0.7188 0.4694 1.531
## Pr(>|z|)
## (Intercept) 9.25e-08 ***
## WHO_participationWHO participation 0.00577 **
## Temperature 0.00230 **
## I(Temperature^2) 0.01310 *
## Prior_self_medicationPatient self medicating 0.01685 *
## WHO_participationWHO participation:Temperature 0.12575
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WHO_pWHOp Tmprtr I(T^2) P__Psm
## WHO_prtWHOp -0.072
## Temperature -0.713 -0.148
## I(Tmprtr^2) 0.597 0.203 -0.956
## Prr_slf_Psm -0.144 0.058 0.007 -0.009
## WHO_pWHOp:T -0.019 -0.858 0.160 -0.243 -0.048Secondary wrong diag
- Wrong diagnosis is additionally significant (p=0.00053).
- Dr. education level with 6 levels did not converge, so we used annotated pooling to combine categories preserving the ordering.
- education level highest (levels 5,6)(p=0.002) and moderate education (levels 3,4) (p=0.04).
fit.t.self.wrd <- glmer(misuse_category2 ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
fit.dred <- glmer(misuse_category2 ~ WHO_participation*Doctor_education_Level+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 2 negative eigenvalues fit.dred <- glmer(misuse_category2 ~ WHO_participation*Doctor_education_Level2+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))Secondary final
- In the multivariate model, education levels were significant (p=0.0002, p=0.006) with an AIC of 821. However age was not (p=0.89) Dropping education feature reduced the AIC to 848 (anova LRT p=0.17).
- Final second model AIC 819.9
final2<- glmer(misuse_category2 ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+age.gc+Doctor_education_Level2+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
summary(final2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## misuse_category2 ~ WHO_participation * Temperature + I(Temperature^2) +
## Prior_self_medication + Wrong_diagnosis + age.gc + Doctor_education_Level2 +
## (1 | Hospital/Doctor.Identifier)
## Data: anti.table
## Control:
## glmerControl(optimizer = "bobyqa", boundary.tol = 3e-05, optCtrl = list(maxfun = 3e+05))
##
## AIC BIC logLik deviance df.resid
## 821.9 878.9 -398.9 797.9 843
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7393 -0.4930 -0.2768 0.1942 3.7759
##
## Random effects:
## Groups Name Variance Std.Dev.
## Doctor.Identifier:Hospital (Intercept) 0.8751 0.9355
## Hospital (Intercept) 0.6014 0.7755
## Number of obs: 855, groups: Doctor.Identifier:Hospital, 134; Hospital, 36
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) -2.02177 0.49012 -4.125
## WHO_participationWHO participation -3.20439 1.20930 -2.650
## Temperature 1.16358 0.39715 2.930
## I(Temperature^2) -0.26006 0.10914 -2.383
## Prior_self_medicationPatient self medicating -0.51603 0.25510 -2.023
## Wrong_diagnosisWrong diagnosis 0.89239 0.24015 3.716
## age.gc -0.01063 0.07873 -0.135
## Doctor_education_Level2Highest education -2.45890 0.66097 -3.720
## Doctor_education_Level2moderate education -0.98107 0.36113 -2.717
## WHO_participationWHO participation:Temperature 0.60733 0.47464 1.280
## Pr(>|z|)
## (Intercept) 3.71e-05 ***
## WHO_participationWHO participation 0.008054 **
## Temperature 0.003392 **
## I(Temperature^2) 0.017182 *
## Prior_self_medicationPatient self medicating 0.043088 *
## Wrong_diagnosisWrong diagnosis 0.000202 ***
## age.gc 0.892618
## Doctor_education_Level2Highest education 0.000199 ***
## Doctor_education_Level2moderate education 0.006594 **
## WHO_participationWHO participation:Temperature 0.200693
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WHO_pWHOp Tmprtr I(T^2) P__Psm Wrn_Wd age.gc D__L2He Dc__L2e
## WHO_prtWHOp -0.049
## Temperature -0.598 -0.155
## I(Tmprtr^2) 0.482 0.210 -0.956
## Prr_slf_Psm -0.127 0.054 0.012 -0.013
## Wrng_dgnsWd -0.385 -0.030 0.031 -0.012 0.063
## age.gc -0.003 0.078 -0.023 0.022 -0.066 -0.015
## Dctr_d_L2He -0.221 0.020 -0.025 0.037 -0.039 -0.102 0.043
## Dctr_dc_L2e -0.432 0.064 -0.021 0.042 -0.036 -0.085 0.026 0.432
## WHO_pWHOp:T -0.015 -0.876 0.162 -0.244 -0.048 -0.007 -0.039 0.018 -0.004final2.noAge<- glmer(misuse_category2 ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+Doctor_education_Level2+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
summary(final2.noAge)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## misuse_category2 ~ WHO_participation * Temperature + I(Temperature^2) +
## Prior_self_medication + Wrong_diagnosis + Doctor_education_Level2 +
## (1 | Hospital/Doctor.Identifier)
## Data: anti.table
## Control:
## glmerControl(optimizer = "bobyqa", boundary.tol = 3e-05, optCtrl = list(maxfun = 3e+05))
##
## AIC BIC logLik deviance df.resid
## 819.9 872.1 -398.9 797.9 844
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7355 -0.4921 -0.2775 0.1937 3.7716
##
## Random effects:
## Groups Name Variance Std.Dev.
## Doctor.Identifier:Hospital (Intercept) 0.8733 0.9345
## Hospital (Intercept) 0.6017 0.7757
## Number of obs: 855, groups: Doctor.Identifier:Hospital, 134; Hospital, 36
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) -2.0221 0.4900 -4.127
## WHO_participationWHO participation -3.1918 1.2044 -2.650
## Temperature 1.1624 0.3970 2.928
## I(Temperature^2) -0.2598 0.1091 -2.381
## Prior_self_medicationPatient self medicating -0.5183 0.2545 -2.036
## Wrong_diagnosisWrong diagnosis 0.8919 0.2401 3.715
## Doctor_education_Level2Highest education -2.4552 0.6601 -3.719
## Doctor_education_Level2moderate education -0.9799 0.3609 -2.715
## WHO_participationWHO participation:Temperature 0.6049 0.4739 1.276
## Pr(>|z|)
## (Intercept) 3.68e-05 ***
## WHO_participationWHO participation 0.008045 **
## Temperature 0.003409 **
## I(Temperature^2) 0.017266 *
## Prior_self_medicationPatient self medicating 0.041710 *
## Wrong_diagnosisWrong diagnosis 0.000204 ***
## Doctor_education_Level2Highest education 0.000200 ***
## Doctor_education_Level2moderate education 0.006631 **
## WHO_participationWHO participation:Temperature 0.201786
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WHO_pWHOp Tmprtr I(T^2) P__Psm Wrn_Wd D__L2He Dc__L2e
## WHO_prtWHOp -0.049
## Temperature -0.598 -0.153
## I(Tmprtr^2) 0.482 0.209 -0.956
## Prr_slf_Psm -0.127 0.060 0.011 -0.011
## Wrng_dgnsWd -0.385 -0.029 0.030 -0.012 0.062
## Dctr_d_L2He -0.221 0.017 -0.024 0.035 -0.037 -0.102
## Dctr_dc_L2e -0.432 0.062 -0.021 0.042 -0.034 -0.085 0.431
## WHO_pWHOp:T -0.015 -0.876 0.161 -0.243 -0.051 -0.007 0.021 -0.003 anova(final2,final2.noAge)Secondary final model
final2<- glmer(misuse_category2 ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+Doctor_education_Level2+(1|Hospital/Doctor.Identifier),
data = anti.table,
family = binomial,
control=glmerControl(optimizer ="bobyqa",boundary.tol=3e-5,optCtrl = list(maxfun=3e5)))
summary(final2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## misuse_category2 ~ WHO_participation * Temperature + I(Temperature^2) +
## Prior_self_medication + Wrong_diagnosis + Doctor_education_Level2 +
## (1 | Hospital/Doctor.Identifier)
## Data: anti.table
## Control:
## glmerControl(optimizer = "bobyqa", boundary.tol = 3e-05, optCtrl = list(maxfun = 3e+05))
##
## AIC BIC logLik deviance df.resid
## 819.9 872.1 -398.9 797.9 844
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7355 -0.4921 -0.2775 0.1937 3.7716
##
## Random effects:
## Groups Name Variance Std.Dev.
## Doctor.Identifier:Hospital (Intercept) 0.8733 0.9345
## Hospital (Intercept) 0.6017 0.7757
## Number of obs: 855, groups: Doctor.Identifier:Hospital, 134; Hospital, 36
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) -2.0221 0.4900 -4.127
## WHO_participationWHO participation -3.1918 1.2044 -2.650
## Temperature 1.1624 0.3970 2.928
## I(Temperature^2) -0.2598 0.1091 -2.381
## Prior_self_medicationPatient self medicating -0.5183 0.2545 -2.036
## Wrong_diagnosisWrong diagnosis 0.8919 0.2401 3.715
## Doctor_education_Level2Highest education -2.4552 0.6601 -3.719
## Doctor_education_Level2moderate education -0.9799 0.3609 -2.715
## WHO_participationWHO participation:Temperature 0.6049 0.4739 1.276
## Pr(>|z|)
## (Intercept) 3.68e-05 ***
## WHO_participationWHO participation 0.008045 **
## Temperature 0.003409 **
## I(Temperature^2) 0.017266 *
## Prior_self_medicationPatient self medicating 0.041710 *
## Wrong_diagnosisWrong diagnosis 0.000204 ***
## Doctor_education_Level2Highest education 0.000200 ***
## Doctor_education_Level2moderate education 0.006631 **
## WHO_participationWHO participation:Temperature 0.201786
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WHO_pWHOp Tmprtr I(T^2) P__Psm Wrn_Wd D__L2He Dc__L2e
## WHO_prtWHOp -0.049
## Temperature -0.598 -0.153
## I(Tmprtr^2) 0.482 0.209 -0.956
## Prr_slf_Psm -0.127 0.060 0.011 -0.011
## Wrng_dgnsWd -0.385 -0.029 0.030 -0.012 0.062
## Dctr_d_L2He -0.221 0.017 -0.024 0.035 -0.037 -0.102
## Dctr_dc_L2e -0.432 0.062 -0.021 0.042 -0.034 -0.085 0.431
## WHO_pWHOp:T -0.015 -0.876 0.161 -0.243 -0.051 -0.007 0.021 -0.003evaluation of secondary model
- The dotplot identified hospital 31, 3 for having a higher association with multiple misuses. Whereas hospitals 10 and 24 had slightly lower than average misuses repeatedly.
myrandef<-ranef(final2,condVar=TRUE)
dotplot(myrandef)## $`Doctor.Identifier:Hospital`
##
## $Hospital
GEE
library("gee");library("geepack");library("nlme")##
## Attaching package: 'nlme'## The following object is masked from 'package:lme4':
##
## lmList## The following object is masked from 'package:dplyr':
##
## collapse library(magrittr);library(dplyr)
anti3 <- anti.table[order(anti.table$Doctor.Identifier), ]
#independence
fit.ind <- gee(misuse_category2 ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+Doctor_education_Level2,
data =anti3 ,
family = binomial(link = "logit"),
corstr = "independence",
id = Doctor.Identifier)## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27## running glm to get initial regression estimate## (Intercept)
## -1.4272419
## WHO_participationWHO participation
## -2.7924947
## Temperature
## 0.7868010
## I(Temperature^2)
## -0.1672340
## Prior_self_medicationPatient self medicating
## -0.6099164
## Wrong_diagnosisWrong diagnosis
## 0.7010512
## Doctor_education_Level2Highest education
## -2.0545555
## Doctor_education_Level2moderate education
## -0.7536521
## WHO_participationWHO participation:Temperature
## 0.5871597gatherEstimates<-function(myFit=NULL){
data<-summary(myFit)$coef
return(as.data.frame(data))
}
ind<-gatherEstimates(myFit=fit.ind)
targs<-c("Estimate","Naive S.E.", "Robust S.E.")
ind<-ind[,targs]
ind$type<-'independence'
ind$corr<-0
##exchangeable.
fit.ex <- gee(misuse_category2 ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+Doctor_education_Level2,
data = anti3,
family =binomial(link = "logit"),
corstr = "exchangeable",
id = Doctor.Identifier)## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
## running glm to get initial regression estimate## (Intercept)
## -1.4272419
## WHO_participationWHO participation
## -2.7924947
## Temperature
## 0.7868010
## I(Temperature^2)
## -0.1672340
## Prior_self_medicationPatient self medicating
## -0.6099164
## Wrong_diagnosisWrong diagnosis
## 0.7010512
## Doctor_education_Level2Highest education
## -2.0545555
## Doctor_education_Level2moderate education
## -0.7536521
## WHO_participationWHO participation:Temperature
## 0.5871597 ex<-gatherEstimates(myFit=fit.ex)
ex<-ex[,targs]
ex$type<-'exchangeable'
ex$corr<-0.1855176
##Autoregressive
fit.ar1<-geeglm(misuse_category2 ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+Doctor_education_Level2,
data = anti3,
family=binomial(link = "logit"),
id=Doctor.Identifier,
corstr='ar1')
ar1<-as.data.frame(summary(fit.ar1)$coef)
ar1[,"Naive S.E."]<-0
ar1[,"Robust S.E."]<-ar1$Std.err
ar1<-ar1[,targs]
ar1$type<-'AR1'
ar1$corr<-0.3057
# unstable
fit.un <- gee(misuse_category2 ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+Doctor_education_Level2,
data = anti3,
family = binomial(link = "logit"),
corstr = "unstructured",
id = Doctor.Identifier)## Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27
## running glm to get initial regression estimate## (Intercept)
## -1.4272419
## WHO_participationWHO participation
## -2.7924947
## Temperature
## 0.7868010
## I(Temperature^2)
## -0.1672340
## Prior_self_medicationPatient self medicating
## -0.6099164
## Wrong_diagnosisWrong diagnosis
## 0.7010512
## Doctor_education_Level2Highest education
## -2.0545555
## Doctor_education_Level2moderate education
## -0.7536521
## WHO_participationWHO participation:Temperature
## 0.5871597un<-gatherEstimates(myFit=fit.un)
un<-un[,targs]
un$type<-'unstructured'
un$corr<-NA
pnorm(coef(summary(fit.un))[,5],lower.tail=TRUE)## (Intercept)
## 9.320087e-06
## WHO_participationWHO participation
## 1.073953e-03
## Temperature
## 9.979258e-01
## I(Temperature^2)
## 1.051462e-02
## Prior_self_medicationPatient self medicating
## 1.279417e-02
## Wrong_diagnosisWrong diagnosis
## 9.997440e-01
## Doctor_education_Level2Highest education
## 1.729813e-05
## Doctor_education_Level2moderate education
## 9.572810e-04
## WHO_participationWHO participation:Temperature
## 8.794587e-01Three level ordinal model
library(VGAM)
anti3<-anti.table
all(anti3$Age==og$age)
all(anti3$Hospital==og$hosp)
all(anti3$Doctor.Identifier==og$doc)
anti3$ordinal.misuse<-1
anti3$ordinal.misuse[which(anti3$Antibiotic_misuse=="misuse of 1 antibiotic")]<-2
anti3$ordinal.misuse[which(anti3$Antibiotic_misuse=="misuse of several antibiotics")]<-3
anti3$ordinal.misuse<-factor(anti3$ordinal.misuse)
library(MASS)
fit.polr<-polr(ordinal.misuse ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+Doctor_education_Level2,
data=anti3,
Hess=TRUE)
summary(fit.polr)
## GEE
library(ordinal)
form<-formula(paste("Antibiotic_misuse ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+Doctor_education_Level2+(1|Hospital)+(1|Doctor.Identifier)",collapse=" "))
fit.re<-clmm(form,
data=anti3)
summary(fit.re)
##relaxed Proportional odds model. (LONGITUDINAL)
# library(mixor)
#fit.mixor.po<-mixor(Antibiotic_misuse ~ WHO_participation*Temperature+I(Temperature**2)+Prior_self_medication+Wrong_diagnosis+Doctor_education_Level2,
# data=anti3,
# id=Doctor.Identifier,
# link='logit',
# nAGQ=20)