######################################################
#Christina Kamis & Jessie West
#Data Expedition
#2/26/21 & 3/3/21
######################################################

## PLEASE SAVE PERMANENT COPY BEFORE EDITING IN TOP RIGHT ##


######################################################
## A couple of notes about R:

# This is a R script. A R script allows us to "save" our work.
# Anything done in the console (below) will disappear once you have
# completed your session in R. Therefore, we want to write all
# code in the R script, so that we can save and re-run our code at anytime.

# You can run a line of code by putting your cursor on that line
# and hitting command + enter. Otherwise, you can highlight a line
# or lines and hit "Run" at the top of the R script. If you do not
# highlight the code before hitting "Run", your entire R script
# will run rather than specific lines of code.

# "#" comments out a line so that the software knows we
# aren't trying to run that line of code, we are simply
# writing a note. This tool is very helpful for annotating our code
# so that we can remember what lines of code mean for future work.

# R is object-oriented, meaning that you must place something
# into an object in order to save that value. To do this
# you simply name the object and then use "<-" to place
# something into that object. For example:

2+2
# The sum of 2+2 is not in an object, so I can't call that value
# back up 

A<-2+2
# I've now placed "2+2" in object called "A"

A
# Now when I call "A" I get the sum of 2+2 ! If you look in the
# "Global Environment" on the right hand side of your screen,
# you will see that "A" is now saved under "Values." I can 
# call "A" anytime I want! 

######################################################
# Before we move forward,
# we must load packages that allow us to give R
# certain commands

#load packages
library(haven)
library(foreign)
library(readr)
library(semPlot)
library(lavaan)
######################################################
# Although the data is in our folder, it has yet to be
# imported into R. To do that we must load the csv.

# load data 
APIMdat <- read_dta("DataExp.dta")
# complete cases (only looking at respondents with no missing data)
APIMdat<-APIMdat[complete.cases(APIMdat),]
# Should have 1,084 Couples

######################################################
## Below are some common R commands for checking out data. You
## won't need these specifically for the assignment, but they can 
## help you better understand the data you are analyzing.


## To view the data set
View(APIMdat)

## I can also get a snapshot of all the variables in my dataset
names(APIMdat)

## To view a specific variable
# Here I am looking at wife cesd at time 1
APIMdat$wifecesdT1
APIMdat$wifecursmokeT2

## To create a table for one variable
table(APIMdat$wifecesdT1)

## To look at frequencies instead of counts
prop.table(table(APIMdat$wifecesdT1))
# You can see from this table that around 55% of wives
# report no depressive symptoms at time 1. 

# Now let's look at husband's depressive symptoms at time 1
prop.table(table(APIMdat$huscesdT1))
# Around 60% of husbands have no depressive symptoms at time 1

## To look at the average value for a variable
mean(APIMdat$wifecesdT1)
mean(APIMdat$huscesdT1)
# The average value is 1.02 for wives and 0.82 for husbands.


######################################################
# If wanting to look at race in your model
# You will need to create new variables 
# for the different race groups rather than using the multiple
# categories variable

table(APIMdat$wiferaracem)
#1= white, 2=Black 3=other

# You can use this code below:
APIMdat$wifewhite<-ifelse(APIMdat$wiferaracem==1,0,1)
APIMdat$huswhite<-ifelse(APIMdat$husraracem==1,0,1)
APIMdat$wifeBlack<-ifelse(APIMdat$wiferaracem==2,0,1)
APIMdat$husBlack<-ifelse(APIMdat$husraracem==1,0,1)

# You should now have four new variables representing identifying
# as "white" or "Black" for husbands and wives.

######################################################
# Here I'm going to create a function that will make the APIM
# process much more user-friendly.
# You just need to run the text below, but don't worry about changing anything

APIMfunc<- function(WifeVar1, WifeVar2, HusVar1, HusVar2) {
  
  
  APIM_1 <- paste("\n\n",
                  HusVar2,  ' ~ a1*',HusVar1,"\n", "
                  ",
                  WifeVar2,' ~ a2*',WifeVar1, "\n","
                  ",
                  HusVar2,' ~ p12*',WifeVar1,"\n","
                  ",
                  WifeVar2,'   ~ p21*',HusVar1,"\n","
                  ",
                  HusVar1 ,' ~ mx1*1',"\n","
                  ",
                  WifeVar1 ,' ~ mx2*1',"\n","
                  ",
                  HusVar2,  ' ~ iy1*1',"\n","
                  ",
                  WifeVar2,' ~ iy2*1',"\n","
                  ",
                  HusVar1, ' ~~ vx1*',HusVar1,"\n","
                  ",
                  WifeVar1, ' ~~ vx2*',WifeVar1,"\n","
                  ",
                  HusVar2,  ' ~~ ve1*',HusVar2,"\n","
                  ",
                  WifeVar2,  ' ~~ ve2*',WifeVar2,"\n","
                  ",
                  HusVar2,' ~~ cy*',WifeVar2,"\n","
                  ",
                  WifeVar1,' ~~ cx*',HusVar1, sep="")
  
  print(APIM_1)
  
 apim1 <- sem(APIM_1,fixed.x=FALSE, data = APIMdat)
 print(summary(apim1))
 
 semPaths(apim1, 
          # general lay-out for a nice APIM:
          fade = F, "est", layout='tree2', rotation = 2, style = "ram",
          intercepts = F, residuals = F, 
          optimizeLatRes = T, curve = 3.1,  
          # labels and their sizes:
          nodeLabels=c(paste(HusVar2), paste(WifeVar2), paste(HusVar1), paste(WifeVar1)),
                       sizeMan=20,  sizeMan2=18,
          # position and size of parameter estimates:
          edge.label.position = 0.45, edge.label.cex=1.5)
}

######################################################
# To get your APIM results, you will need to use the function
# we just created "APIMfunc"

# In order to use this function,
# you will need to write the name of your wife variable from T1,
# wife variable from T2, husband variable from T1, and
# husband variable from T2 inside the () for the APIMfunc

# You will also need to write the name of each variable in quotes
# and separate each quoted variable name with a comma.

# Here is an example:
APIMfunc("wifeadlT1","wifecesdT2","husadlT1","huscesdT2")

# Note that the console will output the "results" from your APIM model
# What you are interested in is under the "Regressions" section of this 
# output.
# Here you will find the "actor" effects (a1 for husband, a2 for wife),
# as well as the "partner" effects (p12 for the effect of wife on husband and
# p21 for the effect of husband on wife).

# We can also look at the "covariances" output to see whether husbands and wives
# are similar at T1 and T2. This would signal health concordance.

# Running APIMfunc will also produce a plot, located at the bottom right hand side
# of R Cloud. 

# Green lines indicate a positive relationship and red lines indicate a negative
# relationship. For this example, a green line indicates that greater disability at T1 
# for a wife is associated with that wife's greater depression at T2. 

# A bold line indicates that this relationship is statistically significant-- in other words,
# we are confident this relationship exists. If the line is thin, that means we
# are not confident in this estimate. For this example, we are confident there is 
# a relationship between husadlT1 and wifecesdT2, but we are NOT confident there is
# a relationship between wifeadlT1 and huscesdT2. 

# Now your turn: 
# REMINDER: you will need to type in the names of variables EXACTLY
# as they appear, in the order presented in the example, and with 
# each variable name in quotes. Otherwise your model will not run/run properly.
APIMfunc()