#peijun

Refer to the questions here.

Q2 Black and White victims

dat <- read.table("http://140.116.183.121/~sheu/lmm/Data/victims.txt", header = T)

# create table to dataframe
tb1 <- with(dat, (table(nvic, race)))
colnames(tb1) <- c("White", "Black")

# run
m0 <- glm.nb(nvic ~ race , data = dat)
summary(m0)

Call:
glm.nb(formula = nvic ~ race, data = dat, init.theta = 0.2023119205, 
    link = log)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.7184  -0.3899  -0.3899  -0.3899   3.5072  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -2.3832     0.1172 -20.335  < 2e-16 ***
race          1.7331     0.2385   7.268 3.66e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(0.2023) family taken to be 1)

    Null deviance: 471.57  on 1307  degrees of freedom
Residual deviance: 412.60  on 1306  degrees of freedom
AIC: 1001.8

Number of Fisher Scoring iterations: 1

              Theta:  0.2023 
          Std. Err.:  0.0409 

 2 x log-likelihood:  -995.7980 
logLik(m0)
'log Lik.' -497.899 (df=3)
AIC(m0)
[1] 1001.798
est <- cbind(Estimate = coef(m0), confint(m0))
Waiting for profiling to be done...
exp(est)
              Estimate      2.5 %   97.5 %
(Intercept) 0.09225413 0.07305984 0.115726
race        5.65841937 3.57782307 9.131638

Q3 borned children

dat <- read.table("http://data.princeton.edu/wws509/datasets/ceb.dat", header = T)
dat <- dat %>% mutate(kids = round(y),os = log(n))

# model
m0 <- glm(formula = kids ~ offset(os), family = "poisson", data = dat)
summary(m0)

Call:
glm(formula = kids ~ offset(os), family = "poisson", data = dat)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-18.6368   -4.4774   -0.8548    2.5292   21.9744  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 1.376346   0.009712   141.7   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 3731.9  on 69  degrees of freedom
Residual deviance: 3731.9  on 69  degrees of freedom
AIC: 4163.3

Number of Fisher Scoring iterations: 5

Q4 Rating gray-scale

dat <- read.table("http://140.116.183.121/~sheu/lmm/Data/grayScale.txt", header =T)
dat$resp <- factor(dat$resp)

# build dataframe
n <- as.data.frame(matrix(dat$count, 12, 5))
names(n) <- c("-2", "-1", "0", "+1", "+2")
rownames(n) <- c("6dark", "5dark", "4dark", "3dark","2dark", "1dark", "1bright", "2bright", "3bright", "4bright", "5bright", "6bright")

# draw
likert(t(n), as.percent = F, main = "", ylab = "Percent")

# model
m0 <- polr(formula = resp ~ contrast, data = dat, weights = count, Hess = T)
summary(m0)
Call:
polr(formula = resp ~ contrast, data = dat, weights = count, 
    Hess = T)

Coefficients:
         Value Std. Error t value
contrast  3.09     0.2354   13.13

Intercepts:
      Value    Std. Error t value 
1|2    -6.1214   0.4958   -12.3470
2|3    -4.5455   0.4275   -10.6325
3|4    -3.2964   0.3617    -9.1136
4|5    -2.4018   0.3300    -7.2793
5|6     0.0089   0.2517     0.0353
6|7     1.0705   0.2700     3.9648
7|8     2.2188   0.3121     7.1088
8|9     3.1776   0.3457     9.1920
9|10    4.3815   0.4068    10.7721
10|11   6.1632   0.4967    12.4088
11|12   6.9703   0.5326    13.0875

Residual Deviance: 594.6589 
AIC: 618.6589

Q6 Domestic Violence

dat <- read.table("http://140.116.183.121/~sheu/lmm/Data/domesticViolence.txt", header = T)
colnames(dat) <- c("Violence", "Advice", "Arrest")

# table
tb1 <- with(dat, table(Violence, Advice))
tb2 <- with(dat, table(Violence, Arrest))

m0 <- glm(Violence ~ Advice + Arrest, family = "binomial", data = dat)
summary(m0)

Call:
glm(formula = Violence ~ Advice + Arrest, family = "binomial", 
    data = dat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.7108  -0.7108  -0.6576  -0.4797   2.1067  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.2470     0.2269  -5.495  3.9e-08 ***
Advice       -0.1744     0.3326  -0.524   0.6001    
Arrest       -0.8571     0.4046  -2.119   0.0341 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 293.66  on 311  degrees of freedom
Residual deviance: 288.59  on 309  degrees of freedom
AIC: 294.59

Number of Fisher Scoring iterations: 4
logLik(m0)
'log Lik.' -144.2949 (df=3)
AIC(m0)
[1] 294.5898
est <- cbind(Estimate = coef(m0), confint(m0))
Waiting for profiling to be done...
exp(est)
             Estimate     2.5 %    97.5 %
(Intercept) 0.2873563 0.1804757 0.4410509
Advice      0.8400000 0.4347063 1.6104854
Arrest      0.4243902 0.1843435 0.9138444
# compare
# Advice
exp(-1.2470 - 0.1744)/(1 + exp(-1.2470 - 0.1744))
[1] 0.1944422
tb1[2,2] / ( tb1[1,2] + tb1[2,2])
[1] 0.1944444
#Arrest
exp(-1.2470 - 0.8571)/(1 + exp(-1.2470 - 0.8571))
[1] 0.108699
tb2[2,2] / ( tb2[1,2] + tb2[2,2])
[1] 0.1086957