R Code

Depression

#summarized frequency data
depression.sum <- data.frame(
  skin_dis = c(1, 1, 0, 0),
  depression = c(1, 0, 1, 0),
  count = c(357, 3176, 58, 1294))
fit1 <- glm(depression ~ skin_dis, weights = count, data = depression.sum,
                family = "binomial")
summary(fit1)
## 
## Call:
## glm(formula = depression ~ skin_dis, family = "binomial", data = depression.sum, 
##     weights = count)
## 
## Deviance Residuals: 
##      1       2       3       4  
##  40.45  -26.01   19.11  -10.65  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.1051     0.1342 -23.135  < 2e-16 ***
## skin_dis      0.9194     0.1454   6.325 2.53e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2840.2  on 3  degrees of freedom
## Residual deviance: 2792.0  on 2  degrees of freedom
## AIC: 2796
## 
## Number of Fisher Scoring iterations: 6
## compare with Wald chi-squared statistic
(coef(fit1) / sqrt(diag(vcov(fit1))))^2 
## (Intercept)    skin_dis 
##   535.22012    40.00574
#crude OR and CI's 
exp(cbind(OR = coef(fit1), confint(fit1)))
## Waiting for profiling to be done...
##                     OR      2.5 %     97.5 %
## (Intercept) 0.04482226 0.03407602 0.05772462
## skin_dis    2.50780639 1.90174378 3.36560836
#comparing to 2x2 table methods from note 5 (not logistic regression)
library(epiR)
## Loading required package: survival
## Package epiR 2.0.52 is loaded
## Type help(epi.about) for summary information
## Type browseVignettes(package = 'epiR') to learn how to use epiR for applied epidemiological analyses
## 
epi.2by2(dat<-c(357, 3176, 58, 1294), method="case.control")
##              Outcome +    Outcome -      Total        Prevalence *        Odds
## Exposed +          357         3176       3533               10.10      0.1124
## Exposed -           58         1294       1352                4.29      0.0448
## Total              415         4470       4885                8.50      0.0928
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Odds ratio                                     2.51 (1.89, 3.33)
## Attrib fraction (est) in the exposed (%)      60.12 (46.81, 70.54)
## Attrib fraction (est) in the population (%)   51.72 (38.44, 62.14)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 42.530 Pr>chi2 = <0.001
## Fisher exact test that OR = 1: Pr>chi2 = <0.001
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units

Anxiety

#summarized frequency data
anxiety.sum <- data.frame(
  skin_dis2 = c(1, 1, 0, 0),
  anxiety = c(1, 0, 1, 0),
  count2 = c(607, 2921, 150, 1200))
fit2 <- glm(anxiety ~ skin_dis2, weights = count2, data = anxiety.sum, 
            family = "binomial")
summary(fit2)
## 
## Call:
## glm(formula = anxiety ~ skin_dis2, family = "binomial", data = anxiety.sum, 
##     weights = count2)
## 
## Deviance Residuals: 
##      1       2       3       4  
##  46.22  -33.21   25.67  -16.81  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.07944    0.08660 -24.013  < 2e-16 ***
## skin_dis2    0.50829    0.09741   5.218 1.81e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 4210.7  on 3  degrees of freedom
## Residual deviance: 4181.4  on 2  degrees of freedom
## AIC: 4185.4
## 
## Number of Fisher Scoring iterations: 5
## compare with Wald chi-squared statistic
(coef(fit2) / sqrt(diag(vcov(fit2))))^2
## (Intercept)   skin_dis2 
##   576.62648    27.22792
#crude OR and CI's 
exp(cbind(OR = coef(fit2), confint(fit2)))
## Waiting for profiling to be done...
##                   OR     2.5 %    97.5 %
## (Intercept) 0.125000 0.1050788 0.1475891
## skin_dis2   1.662444 1.3771519 2.0180425
#comparing to 2x2 table methods from note 5 (not logistic regression)
epi.2by2(dat<-c(607, 2921, 150, 1200), method="case.control")
##              Outcome +    Outcome -      Total        Prevalence *        Odds
## Exposed +          607         2921       3528                17.2       0.208
## Exposed -          150         1200       1350                11.1       0.125
## Total              757         4121       4878                15.5       0.184
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Odds ratio                                     1.66 (1.37, 2.01)
## Attrib fraction (est) in the exposed (%)      39.84 (27.01, 50.65)
## Attrib fraction (est) in the population (%)   31.95 (20.86, 41.49)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 27.658 Pr>chi2 = <0.001
## Fisher exact test that OR = 1: Pr>chi2 = <0.001
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units

Suicidal Ideation

#summarized frequency data
suic.sum <- data.frame(
  skin_dis3 = c(1, 1, 0, 0),
  suicide = c(1, 0, 1, 0),
  count3 = c(451, 3184, 88, 1271))
fit3 <- glm(suicide ~ skin_dis3, weights = count3, data = suic.sum,
                family = "binomial")
summary(fit3)
## 
## Call:
## glm(formula = suicide ~ skin_dis3, family = "binomial", data = suic.sum, 
##     weights = count3)
## 
## Deviance Residuals: 
##      1       2       3       4  
##  43.39  -29.04   21.95  -13.04  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.6702     0.1102 -24.224  < 2e-16 ***
## skin_dis3     0.7158     0.1212   5.907 3.47e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3417.5  on 3  degrees of freedom
## Residual deviance: 3377.9  on 2  degrees of freedom
## AIC: 3381.9
## 
## Number of Fisher Scoring iterations: 6
## compare with Wald chi-squared statistic
(coef(fit3) / sqrt(diag(vcov(fit3))))^2 
## (Intercept)   skin_dis3 
##   586.81929    34.89801
#crude OR and CI's 
exp(cbind(OR = coef(fit3), confint(fit3)))
## Waiting for profiling to be done...
##                     OR      2.5 %     97.5 %
## (Intercept) 0.06923682 0.05539221 0.08537513
## skin_dis3   2.04581501 1.62164811 2.60900446
#comparing to 2x2 table methods from note 5 (not logistic regression)
epi.2by2(dat<-c(451, 3184, 88, 1271), method="case.control")
##              Outcome +    Outcome -      Total        Prevalence *        Odds
## Exposed +          451         3184       3635               12.41      0.1416
## Exposed -           88         1271       1359                6.48      0.0692
## Total              539         4455       4994               10.79      0.1210
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Odds ratio                                     2.05 (1.61, 2.59)
## Attrib fraction (est) in the exposed (%)      51.11 (37.83, 61.90)
## Attrib fraction (est) in the population (%)   42.77 (30.33, 52.99)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 36.150 Pr>chi2 = <0.001
## Fisher exact test that OR = 1: Pr>chi2 = <0.001
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units