hw2

1 data input

dta2 <- read.table("cantabrian.txt", h=T)
names(dta2) <- c("ID","Health","Gender","Type")
head(dta2)

##   ID Health Gender Type
## 1  1      0      M   GP
## 2  1      0      M  GHQ
## 3  2      0      M   GP
## 4  2      0      M  GHQ
## 5  3      0      M   GP
## 6  3      0      M  GHQ

2 data management

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, union

dta2t <- dta2 %>% 
  mutate(Health = factor(Health, 
                         levels=0:1, 
                         labels=c("Illness","Otherwise")),
         Gender = factor(Gender, levels = c("M", "F")),
         Type = factor(Type, 
                       levels=c("GP","GHQ")))
str(dta2t)

## 'data.frame':    1646 obs. of  4 variables:
##  $ ID    : int  1 1 2 2 3 3 4 4 5 5 ...
##  $ Health: Factor w/ 2 levels "Illness","Otherwise": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Gender: Factor w/ 2 levels "M","F": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Type  : Factor w/ 2 levels "GP","GHQ": 1 2 1 2 1 2 1 2 1 2 ...

str(dta2)

## 'data.frame':    1646 obs. of  4 variables:
##  $ ID    : int  1 1 2 2 3 3 4 4 5 5 ...
##  $ Health: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Gender: chr  "M" "M" "M" "M" ...
##  $ Type  : chr  "GP" "GHQ" "GP" "GHQ" ...

3 Descriptive

ftable(dta2t, row.vars = c('Gender','Health'), col.vars = c('Type'))

##                  Type  GP GHQ
## Gender Health                
## M      Illness        287 244
##        Otherwise       37  80
## F      Illness        420 306
##        Otherwise       79 193

prop.table(dta2p <- (ftable(dta2t, row.vars = c('Gender','Health'), col.vars = c('Type'))), 2)

##                  Type         GP        GHQ
## Gender Health                              
## M      Illness        0.34872418 0.29647631
##        Otherwise      0.04495747 0.09720535
## F      Illness        0.51032807 0.37181045
##        Otherwise      0.09599028 0.23450790

女性精神病的盛行率高於男性

4 GEE

library(geepack)
summary(m0 <- geeglm(Health ~ Gender + Type + Gender*Type, data=dta2, id = ID, corstr="exchangeable", family = binomial))

## 
## Call:
## geeglm(formula = Health ~ Gender + Type + Gender * Type, family = binomial, 
##     data = dta2, id = ID, corstr = "exchangeable")
## 
##  Coefficients:
##                Estimate  Std.err   Wald Pr(>|W|)    
## (Intercept)    -0.46089  0.09192 25.141 5.33e-07 ***
## GenderM        -0.65425  0.15826 17.089 3.57e-05 ***
## TypeGP         -1.20991  0.12651 91.462  < 2e-16 ***
## GenderM:TypeGP  0.27649  0.22828  1.467    0.226    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = exchangeable 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)        1 0.07865
##   Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha   0.2975 0.04744
## Number of clusters:   823  Maximum cluster size: 2

sjPlot::tab_model(m0, show.obs=F, show.ngroups = F)

	Health
Predictors	Odds Ratios	CI	p
(Intercept)	0.63	0.53 – 0.76	<0.001
Gender [M]	0.52	0.38 – 0.71	<0.001
Type [GP]	0.30	0.23 – 0.38	<0.001
Gender [M] * Type [GP]	1.32	0.84 – 2.06	0.226

女性患病率是男性的0.52倍

GHQ的測出具精神病是GP的0.3倍

性別與檢測方式(GP or GHQ)沒有交互作用