#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
#@@@ logistic regression in R @@@
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

# To assess how covariates affect assortative mating.
# independent variables:
# am or pm- a varaible. character data, essentially binary.
# male holding number- a variable. categorical, essentially binary.
# population- 1 or 2. categorical, essentially binary.
# wing clip effect. 

# Dependent variable:
# frequency of assortative mating: the ratio of assort/total matings. numeric continuous data.

#Thus for a logistic regression:
# predictor/independent variables = all covariates
# response/dependent variable = counts of assort and disassort per experiment


#step 1: upload data in r:
setwd("C:/Program Files/RStudio/Man R-directory")

all.data = read.csv("logsitic regression data for R analysis on 5.15.15.csv")

#step 2: manipulate all.data to obtian counts of assort
# and disassort matings and frequencies of each.the include next
# variable of wing clip effect.
summary(all.data) 
str(all.data)
all.data$a.freq <- (all.data$assort/all.data$N)
all.data$d.freq <- ((all.data$N-all.data$assort)/all.data$N)
all.data$disassort <- (all.data$N-all.data$assort)

all.data$wingclipeffect <- c("no")
all.data[16:19,11] <- "yes"

# here turn on library plyr
library(plyr)
library(reshape)
library(reshape2)

all.data <- all.data[,-(8:9)]
all.data[-21,]
all.data <- all.data[-21,]
melt.data.frame(all.data, id.var = c("date","strain","male.hold",
                                     "population","time.day","N","wingclipeffect"))
all.data <- melt.data.frame(all.data, id.var = c("date","strain","male.hold",
                                                 "population","time.day","N","wingclipeffect"))
all.data$freq <- (all.data$value/all.data$N)
mode(all.data$freq)

#step 3: split all.data into two seperate dataframes based on genotype.
order(all.data$strain, all.data$variable, all.data$date)
index <- order(all.data$strain, all.data$variable, all.data$date)
freq.data <-all.data[index,]

all.ral <- as.data.frame(freq.data[1:16,])
can.s <- as.data.frame(freq.data[17:40,])

# step 4: run a regression. 
# All ral
allral.regression <- glm(variable~value + time.day + population,
                      data = all.ral, family = binomial)
summary(allral.regression)


cans.regression <- glm(variable~value + time.day + population
                       + male.hold + wingclipeffect,data = can.s, family = binomial)
summary(cans.regression)

#Step 5: confidence intervals
confint(allral.regression)
confint(cans.regression)