DS_assignment7

C ode Setup.
Load required packages.

library(tidyverse)
library(keras)
library(tensorflow)

Question 1:

Weighted least squares. When input vector x is weighted by a vector w , the least squares solution is obtained as :
\[ \mu_{optimal } = \frac{\sum_{i=1}^{n} (w_{i} *x_{i})} {\sum_{i=1}^{n} w_{i}} \]

# Question 1
x = c(1.57, 1.25, 2.80, 0.43)
w =c(2, 2, 1, 1)
optimal_mu = sum(x*w)/ sum(w)
print(paste("Optimal value ", optimal_mu))

## [1] "Optimal value  1.47833333333333"

Question 2: Linear Regression

# Question 2 
# Load dataset
shhs =  read.delim("/Users/archanabalan/Documents/Spring Semester/DS/shhs.txt")
# Add a new column as: log(rdi4p + 1)
log_add = function(x){
  log(x + 1)
}
shhs$log_rdi4p = unlist(lapply(shhs$rdi4p, log_add))
# Fit a linear regression modl
model_q2 = lm(log_rdi4p ~ waist + age_s1 + gender + bmi_s1, data = shhs,na.action=na.exclude)
model_coef = coef(model_q2)
print("The model coefficients are:" )

## [1] "The model coefficients are:"

print(model_coef)

##  (Intercept)        waist       age_s1       gender       bmi_s1 
## -2.281391771  0.007058175  0.019982473  0.517820827  0.063067594

Question 3: Predicting using linear model

Inputs: waist =100, age = 60, gender = 0 (male) , bmi = 30

# Question 3
# waist =100, age = 60, gender = 0 (male) , bmi = 30
input_data = c(100,60,0,30)
names(input_data)= c("waist", "age_s1", "gender", "bmi_s1")
out_log_rdi4p = predict(model_q2, newdata = as.data.frame(t(input_data)))
out_rdi4p = exp(out_log_rdi4p)-1
print(paste("The estimated value of log(rdi4p + 1) is ", out_log_rdi4p))

## [1] "The estimated value of log(rdi4p + 1) is  1.51540196188145"

print(paste("Thus, the  estimated value of rdi4p  is ", out_rdi4p))

## [1] "Thus, the  estimated value of rdi4p  is  3.55125018830316"

Question 4: Interpretation of model coefficients

Intercept: It is the value of log(rdi4p + 1) when all other values are zero.
Waist: When values of all other input variables are fixed the value of log(rdi4p + 1) would increase by 0.007058175 for a unit increase in waist.
Age:When values of all other input variables are fixed the value of log(rdi4p + 1) would increase by 0.019982473 for a unit increase in age.
Gender:When values of all other input variables are fixed the value of log(rdi4p + 1) would increase by 0.517820827 for a unit increase in age.This implies that if the person is male(gender= 0) there will be no change in output but if the person is feamle(gender = 1) the ouput will increase by 0.517820827 if all other input variables remain constant.
BMI:When values of all other input variables are fixed the value of log(rdi4p + 1) would increase by 0.063067594 for a unit increase in BMI.

Question 5: Logistic Regression

# Question 5  
# Add a new column to the dataframe defining the categories of rdi4p , based on the definitions provided
# Categorical Variables
# function to find the category of rdi4p
category_rdi4p = function(rdi4p){
  if (rdi4p< 7){ category = "Normal Range"}
  else if (rdi4p >= 7 & rdi4p <15 ){ category = "Mild Sleep Apnea "}
  else if (rdi4p >= 15 & rdi4p <30){ category = " Sleep Apnea "}
  else  if (rdi4p >= 30){ category = " Severe Sleep Apnea "}
}
# Computing the categories
shhs$rdi4p_categories = unlist(lapply(shhs$rdi4p, category_rdi4p))

# Logistic regression
model_q5 = glm(HTNDerv_s1~ rdi4p_categories, data = shhs, family = "binomial")
model_q5_coef = coef(model_q5)
summary(model_q5)

## 
## Call:
## glm(formula = HTNDerv_s1 ~ rdi4p_categories, family = "binomial", 
##     data = shhs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3328  -0.9805  -0.9805   1.2007   1.3880  
## 
## Coefficients:
##                                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                         0.3581     0.1086   3.296 0.000979 ***
## rdi4p_categories Sleep Apnea       -0.3703     0.1338  -2.767 0.005665 ** 
## rdi4p_categoriesMild Sleep Apnea   -0.4126     0.1242  -3.321 0.000897 ***
## rdi4p_categoriesNormal Range       -0.8407     0.1138  -7.390 1.47e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 7921.7  on 5803  degrees of freedom
## Residual deviance: 7824.8  on 5800  degrees of freedom
## AIC: 7832.8
## 
## Number of Fisher Scoring iterations: 4

print("The model coefficients are:" )

## [1] "The model coefficients are:"

print(model_q5_coef)

##                       (Intercept)     rdi4p_categories Sleep Apnea  
##                         0.3580629                        -0.3702954 
## rdi4p_categoriesMild Sleep Apnea       rdi4p_categoriesNormal Range 
##                        -0.4126219                        -0.8407338

Question 6: Logistic Regression with additional variables

# Question 6  

model_q6 = glm(HTNDerv_s1~ rdi4p_categories + waist + age_s1 + gender + bmi_s1+ smokstat_s1, data = shhs, family = "binomial")
model_q6_coef = coef(model_q6)
summary(model_q6)

## 
## Call:
## glm(formula = HTNDerv_s1 ~ rdi4p_categories + waist + age_s1 + 
##     gender + bmi_s1 + smokstat_s1, family = "binomial", data = shhs)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9029  -1.0234  -0.6827   1.1240   2.1408  
## 
## Coefficients:
##                                    Estimate Std. Error z value Pr(>|z|)
## (Intercept)                       -5.122034   0.338140 -15.148  < 2e-16
## rdi4p_categories Sleep Apnea      -0.341663   0.143284  -2.385 0.017102
## rdi4p_categoriesMild Sleep Apnea  -0.303458   0.133779  -2.268 0.023308
## rdi4p_categoriesNormal Range      -0.452920   0.126855  -3.570 0.000356
## waist                              0.004466   0.003620   1.234 0.217298
## age_s1                             0.055657   0.002937  18.952  < 2e-16
## gender                             0.012453   0.065102   0.191 0.848310
## bmi_s1                             0.044391   0.009553   4.647 3.37e-06
## smokstat_s1                       -0.023995   0.031149  -0.770 0.441097
##                                      
## (Intercept)                       ***
## rdi4p_categories Sleep Apnea      *  
## rdi4p_categoriesMild Sleep Apnea  *  
## rdi4p_categoriesNormal Range      ***
## waist                                
## age_s1                            ***
## gender                               
## bmi_s1                            ***
## smokstat_s1                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 7380.2  on 5383  degrees of freedom
## Residual deviance: 6837.8  on 5375  degrees of freedom
##   (420 observations deleted due to missingness)
## AIC: 6855.8
## 
## Number of Fisher Scoring iterations: 4

print("The model coefficients are:" )

## [1] "The model coefficients are:"

print(model_q6_coef)

##                       (Intercept)     rdi4p_categories Sleep Apnea  
##                      -5.122033742                      -0.341662580 
## rdi4p_categoriesMild Sleep Apnea       rdi4p_categoriesNormal Range 
##                      -0.303457620                      -0.452920307 
##                             waist                            age_s1 
##                       0.004465803                       0.055657331 
##                            gender                            bmi_s1 
##                       0.012452504                       0.044390888 
##                       smokstat_s1 
##                      -0.023995221

Results interpretation: The model chooses only 3 categorical variables (Sleep Apnea, Mild Sleep Apnea, Normal Range). Thus the presence or absence of Severe Sleep Apnea does not have any impact on hypertension.

Intercept: If all input variables are zero the value of hypertension is -5.122033742.

When all other input variables are kept constant a unit increase in Sleep Apnea relatively decreases hypertension by 0.341662580.
When all other input variables are kept constant a unit increase in Mild Sleep Apnea relatively decreases hypertension by 303457620.
When all other input variables are kept constant a unit increase in Normal Range relatively decreases hypertension by 0.452920307.
When all other input variables are kept constant a unit increase in waist relatively increases hypertension by 0.004465803.
When all other input variables are kept constant a unit increase in age relatively increases hypertension by 0.055657331.
When all other input variables are kept constant a unit increase in smokstat_s1 relatively decreases hypertension by 0.023995221.

Question 7: Repeating linear regression model with neural network

Estimation of log(rdi4p+1) using regression variables age, waist, gender and bmi (as given in question 2)

# Question 7 
set.seed(123)
# Choose the variables of interest
data_nn = shhs %>% select(log_rdi4p ,waist , age_s1 ,gender , bmi_s1)
data_nn = na.omit(data_nn)
# Define output and input variables
data_nn_y = data_nn$log_rdi4p
data_nn_x = data_nn %>% select(-log_rdi4p) %>% as.matrix()
# Split training (70% ) and testing (30%)
train_index = sample(c(TRUE, FALSE), length(data_nn_y), replace = TRUE, prob = c(.7, .3))
train_y = data_nn_y[train_index]
test_y = data_nn_y[!train_index]

train_x = data_nn_x[train_index, ] 
test_x = data_nn_x[!train_index, ] 
# Scale and centre the input data
mean    <- colMeans(train_x)
std     <- apply(train_x, 2, sd)
train_x <- scale(train_x, center = mean, scale = std)
test_x <- scale(test_x, center = mean, scale = std)

# Set up the model
model <- keras_model_sequential() %>%
  # network architecture
  layer_dense(units = 16,activation = "relu",use_bias = TRUE,  input_shape = ncol(train_x)) %>%
  layer_dense(units = 4,activation = "relu") %>%
  layer_dense(units = 2,activation = "relu") %>%
  layer_dense(units = 1)
# Print model summary
summary(model)

## ___________________________________________________________________________
## Layer (type)                     Output Shape                  Param #     
## ===========================================================================
## dense (Dense)                    (None, 16)                    80          
## ___________________________________________________________________________
## dense_1 (Dense)                  (None, 4)                     68          
## ___________________________________________________________________________
## dense_2 (Dense)                  (None, 2)                     10          
## ___________________________________________________________________________
## dense_3 (Dense)                  (None, 1)                     3           
## ===========================================================================
## Total params: 161
## Trainable params: 161
## Non-trainable params: 0
## ___________________________________________________________________________

# backpropagation
model %>% compile(
    optimizer = optimizer_rmsprop(),
    loss = "mse",
    metrics =  list("mean_absolute_error")
    )

# Training the model
learn <- model %>% fit(
  x = train_x,
  y = train_y,
  epochs = 40,
  batch_size = 32,
  validation_split = .2,
  verbose = 1

)
# Print training output 
print("Trained Model: ")

## [1] "Trained Model: "

learn

## Trained on 3,043 samples, validated on 761 samples (batch_size=32, epochs=40)
## Final epoch (plot to see history):
##                val_loss: 0.7403
## val_mean_absolute_error: 0.6979
##                    loss: 0.8391
##     mean_absolute_error: 0.7453

## compare with the test data
pred_nn = model %>% predict(test_x)

## compare with ordinary regression
fit = lm(train_y ~ train_x)
pred_lm =  cbind(1, test_x) %*% coef(fit)
cor_nn = cor(pred_nn, test_y)
cor_lm = cor(pred_lm, test_y)
# 
print(paste("Correlation with the test set for neural network, ",cor_nn))

## [1] "Correlation with the test set for neural network,  0.518543024505331"

print(paste("Correlation with the test set for linear regression, ",cor_lm))

## [1] "Correlation with the test set for linear regression,  0.514075641770967"

**Note: Reported is the best performing neural network after tuning for hyperparamters such as epochs and multiple layers

Question 8: Repeating logistic regression model with neural network

Estimation of HTNDerv_s1 using regression variables rdi4p (categorical), age, waist, gender , bmi, smoke_stat_s1 (as given in question 6)

# Question 8
set.seed(123)
# Choose the variables of interest
data_nn = shhs %>% select(rdi4p, waist, HTNDerv_s1, gender, age_s1, bmi_s1)
data_nn = na.omit(data_nn)
# Convert output to a binary output variable
data_nn_y = data_nn$HTNDerv_s1
data_nn_y_binary  = to_categorical(data_nn_y,2)
 # Define  input variables
data_nn_x = data_nn %>% select(-HTNDerv_s1) %>% as.matrix()
# Split training (70% ) and testing (30%)
train_index = sample(c(TRUE, FALSE), dim(data_nn_x )[1], replace = TRUE, prob = c(.7, .3))

train_y = data_nn_y_binary [train_index, ]
train_x = data_nn_x[train_index, ] %>% scale()

test_y = data_nn_y_binary [!train_index, ]
test_x = data_nn_x[!train_index, ] %>% scale()
# Scale and centre the input data
mean    <- colMeans(train_x)
std     <- apply(train_x, 2, sd)
train_x <- scale(train_x, center = mean, scale = std)
test_x <- scale(test_x, center = mean, scale = std)

# Set up the model
model = keras_model_sequential() %>%
  layer_dense(units = 2^8, activation = "relu", input_shape = ncol(train_x)) %>% 
  layer_dropout(rate = 0.4) %>%
  layer_dense(units = 2^4, activation = "relu") %>%
  layer_dropout(rate = 0.3) %>%
  layer_dense(units = 2, activation = "softmax")
summary(model)

## ___________________________________________________________________________
## Layer (type)                     Output Shape                  Param #     
## ===========================================================================
## dense_4 (Dense)                  (None, 256)                   1536        
## ___________________________________________________________________________
## dropout (Dropout)                (None, 256)                   0           
## ___________________________________________________________________________
## dense_5 (Dense)                  (None, 16)                    4112        
## ___________________________________________________________________________
## dropout_1 (Dropout)              (None, 16)                    0           
## ___________________________________________________________________________
## dense_6 (Dense)                  (None, 2)                     34          
## ===========================================================================
## Total params: 5,682
## Trainable params: 5,682
## Non-trainable params: 0
## ___________________________________________________________________________

# backpropagation
model %>% compile(
  optimizer = optimizer_rmsprop(),
  loss = "categorical_crossentropy",
  metrics =  list("accuracy")
)
# Training the model
learn <- model %>% fit(
  x = train_x,
  y = train_y,
  epochs = 40,
  batch_size = 32,
  validation_split = .2,
  verbose = 1
)

# Print the trained model
learn

## Trained on 3,043 samples, validated on 761 samples (batch_size=32, epochs=40)
## Final epoch (plot to see history):
## val_loss: 0.6367
##  val_acc: 0.6216
##     loss: 0.6281
##      acc: 0.6477

## compare with the test data
print(" Model performance on test data")

## [1] " Model performance on test data"

acc_nn=model %>% evaluate(test_x, test_y)

## compare with ordinary logistic regression
fit = glm(HTNDerv_s1 ~ ., family = "binomial", data = data_nn, subset = train_index)
pred_lg =  (predict(fit, data_nn[!train_index,], type = "response") > 0.5) * 1
 
ptab2 = table(pred_lg, test_y[,1])
acc_lg = sum(diag(ptab2)) / sum(ptab2)

print(paste("Correlation with the test set for neural network, ",acc_nn$acc))

## [1] "Correlation with the test set for neural network,  0.609801488981649"

print(paste("Correlation with the test set for logistic regression, ",acc_lg))

## [1] "Correlation with the test set for logistic regression,  0.387096774193548"