#####################################################################
#
# Construction and analysis of Table 11.2 in Agresti, for pp 506-507
#
#####################################################################

Freq <- c(16, 13,  9,  3, 14,  4, 15,  6,
          31,  0,  6,  0, 22,  2,  9,  0,
           2,  2,  8,  9,  9, 15, 27, 28,
           7,  2,  5,  2, 31,  5, 32,  6)

Diag <- c(rep("Mild",16),rep("Severe",16))

Tx <- rep(c(rep("Std",8),rep("New",8)),2)

R1 <- rep(c(1,1,1,1,0,0,0,0),4)  # 1 = normal, 0 = abnormal

R2 <- rep(c(1,1,0,0,1,1,0,0),4)

R3 <- rep(c(1,0,1,0,1,0,1,0),4)

table11.2.short <- data.frame(Diag,Tx,R1,R2,R3,Freq)

#
# table11.2.long:
#
# case = patient (case) identifier
#    Diag = severity of depression (diagnosis)
#    Drug = drug tx (std or new)
#    Time = time (1, 2, or 4 weeks) [though 0,1,2 is another possibility]
#    y = outcome (1=normal,0=abnormal)
#
table11.2.long <- data.frame(case=c(0,0), Diag=c("Mild","Severe"),
                             Drug=c("Std","New"), Time=0, y=0)
n <- 0
for (i in 1:32) {
   if (Freq[i]>0) for (k in 1:Freq[i]) {
      n <- n + 1
      for (j in 1:3) {
         table11.2.long <- rbind(table11.2.long,
	   c(case=n, Diag=Diag[i], Drug=Tx[i],
########## Time=2^(j-1),y=c(R1[i],R2[i],R3[i])[j]))
           Time=j,y=c(R1[i],R2[i],R3[i])[j]))
      }
   }
}
table11.2.long <- table11.2.long[-(1:2),]
rownames(table11.2.long) <- paste(1:dim(table11.2.long)[1])

table11.2.long$y <- as.numeric(table11.2.long$y)
table11.2.long$Time <- as.numeric(table11.2.long$Time)
table11.2.long$Diag <- ifelse(table11.2.long$Diag=="Severe",1,0)
table11.2.long$Drug <- ifelse(table11.2.long$Drug=="New",1,0)
table11.2.long$case <- as.factor(table11.2.long$case)

rm(Freq,Diag,Tx,R1,R2,R3)

###################################################################

# Fixed effects model - the usual logistic regression model

fixed.fit <- glm(y ~ Diag + Drug*Time,data=table11.2.long,family=binomial)

summary(fixed.fit)

# Call:
# glm(formula = y ~ Diag + Drug * Time, family = binomial, data = table11.2.long)
# 
# Deviance Residuals: 
#     Min       1Q   Median       3Q      Max  
# -2.4352  -1.0220   0.3254   0.9915   1.8009  
# 
# Coefficients:
#             Estimate Std. Error z value Pr(>|z|)    
# (Intercept)  -0.5104     0.2561  -1.993  0.04627 *  
# Diag         -1.3139     0.1464  -8.974  < 2e-16 ***
# Drug         -1.0770     0.3847  -2.800  0.00512 ** 
# Time          0.4824     0.1148   4.204 2.63e-05 ***
# Drug:Time     1.0174     0.1888   5.389 7.08e-08 ***
# ---
# Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 
#
# (Dispersion parameter for binomial family taken to be 1)
# 
#     Null deviance: 1411.9  on 1019  degrees of freedom
# Residual deviance: 1161.9  on 1015  degrees of freedom
# AIC: 1171.9
# 
# Number of Fisher Scoring iterations: 4
#


###################################################################

# Mixed effects model: the logistic-normal model with a random
# intercept for each person (case) in the file

library(glmmML)
rand.incpt <- glmmML(y ~ Diag + Drug*Time,cluster=case,
   data=table11.2.long,family=binomial)
rand.incpt

# Call:  glmmML(formula = y ~ Diag + Drug * Time, 
#        data = table11.2.long,cluster = case, family = binomial) 
#
#
#                coef se(coef)      z Pr(>|z|)
# (Intercept) -0.5108   0.2567 -1.990 4.66e-02
# Diag        -1.3152   0.1546 -8.505 0.00e+00
# Drug        -1.0781   0.3871 -2.785 5.35e-03
# Time         0.4828   0.1160  4.164 3.13e-05
# Drug:Time    1.0184   0.1924  5.293 1.20e-07
#
# Standard deviation in mixing distribution:  0.06581 
# Std. Error:                                 1.242 
#
# Residual deviance: 1162 on 1014 degrees of freedom      AIC: 1174 
# 

##############################################################

library(MASS)
 
rand.incptPQL <- glmmPQL(y ~ Diag + Drug*Time, 
   random = ~ 1 | case, family=binomial, data=table11.2.long)
 
summary(rand.incptPQL)
 
# Linear mixed-effects model fit by maximum likelihood
#  Data: table11.2.long 
#        AIC      BIC    logLik
#   4611.546 4646.039 -2298.773
# 
# Random effects:
#  Formula: ~1 | case
#         (Intercept)  Residual
# StdDev:   0.3151588 0.9719163
# 
# Variance function:
#  Structure: fixed weights
#  Formula: ~invwt 
# Fixed effects: y ~ Diag + Drug * Time 
#                  Value Std.Error  DF   t-value p-value
# (Intercept) -0.5089348 0.2515715 678 -2.023022  0.0435
# Diag        -1.3170725 0.1469893 337 -8.960329  0.0000
# Drug        -1.0824002 0.3770374 337 -2.870803  0.0044
# Time         0.4825385 0.1118404 678  4.314528  0.0000
# Drug:Time    1.0209475 0.1843255 678  5.538828  0.0000
#  Correlation: 
#           (Intr) Diag   Drug   Time  
# Diag      -0.227                     
# Drug      -0.641  0.035              
# Time      -0.876 -0.090  0.595       
# Drug:Time  0.562 -0.080 -0.921 -0.595
# 
# Standardized Within-Group Residuals:
#        Min         Q1        Med         Q3        Max 
# -4.1651301 -0.8139819  0.2299254  0.8120307  2.1346920 
# 
# Number of Observations: 1020
# Number of Groups: 340 
# 

# Let's look at the fitted odds ratios...

multicenter <- cbind(multicenter,fitted=
  unique(predict(glmmPQL.fit2,
  multicenter.big,type="response")))




multicenter <- multicenter[,-6]