---
title: "350.wk6.programming"
author: "William Revelle"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
options(width=110) 
```

```{r}
#first make psych and psych tools active
library(psych)
library(psychTools)
```


# Modifying (Improving) a function

Some functions are derived from basic principles.  Other functions are merely modifications (improvements - extensions) of pre existing functions.

## Consider `setCor.diagram`.

The `lmCor` function mimics `lm` for doing regressions.  `lmCor` works from correlation matrices as well.  `lmCor` also draws a figure (using `setCor.diagram') showing the path models.  `lm` does not do this.  Can we modify `setCor.diagram' to do this?

## First compare the output of `lm` and `lmCor`.  

We use the attitude data set built into R.

From the help page:

From a survey of the clerical employees of a large financial organization, the data are aggregated from the questionnaires of the approximately 35 employees for each of 30 (randomly selected) departments. The numbers give the percent proportion of favourable responses to seven questions in each department.

```{r}
describe(attitude)

mod.lm <- lm(rating ~ complaints + privileges + raises, data=attitude)
summary(mod.lm)
```

Then do the same using `lmCor`.  Although written originally to just work from correlation matrices, `lmCor` was extended to work from the raw data.

```{r}
mod.1 <- lmCor(rating ~ complaints + privileges + raises, data=attitude)
mod.1
```

Why are these different?  `lmCor` by default standardizes the variables before doing the analysis. We can turn this off by saying std=FALSE.

```{r}
mod.2 <- lmCor(rating ~ complaints + privileges + raises, data=attitude,std=FALSE)
mod.2
```

## examine the names and classes of mod2 and mod.lm

```{r}
class(mod.lm)
class(mod.2)
names(mod.lm)
names(mod.2)
```

### compare some of these elements

```{r}
mod.lm$coefficients
mod.2$coefficients
```


# `lmCor` calls setCor.diagram

Lets look at that function to see how it works and figure out how to use it to do diagram from the `lm` output.

```{r}
setCor.diagram
```
Lets examine some of the parameters being passed to the function

```{r}
mod.2$x.matrix
mod.2$y.matrix
mod.2$resid
```

To make these numbers slightly more understandable, lets use the standardized regression (mod.1).

```{r}
mod.1$x.matrix
mod.1$y.matrix
mod.1$resid
```

## A revised setCor.diagram function

We have to choose whether we are doing `lmCor` or `lm` analyses.

Then, if `lm` we need to find the covariance matrix from the data, rather than have it passed to us.


```{r}
lmDiagram <- function(sc,main="Regression model",digits=2,show=FALSE,cex=1,l.cex=1,...) { 
if(missing(l.cex)) l.cex <- cex  

beta <- round(sc$coefficients,digits)
#two cases
#the normal case, the model comes from lmCor
#the 2nd case, the model comes from lm
if(length(class(sc)) ==1) {lm<- TRUE} else { lm <- FALSE}
if(!lm) {#case 1
if(rownames(beta)[1] %in% c("(Intercept)", "intercept*")) {intercept <- TRUE} else {intercept <- FALSE}
   
x.matrix <- round(sc$x.matrix,digits)

y.matrix <- round(sc$y.matrix,digits)
y.resid <- round(sc$resid,digits)
x.names <- rownames(beta)

if(intercept){ x.matrix <- x.matrix[-intercept,-intercept]
 x.names <- x.names[-intercept] 
 }
beta <- beta[-intercept,,drop=FALSE]
y.names <- colnames(beta)
} else { x.matrix <- round(cov(sc$model[-1]), digits=digits)  #case 2  from lm
       x.names <- colnames(sc$model)[-1]
       beta <- as.matrix(beta[-1]) #make the vector into a matrix, drop the intercept
      y.names <-  colnames(sc$model)[1]
      }
      
nx <- length(x.names)
ny <- length(y.names) 
top <- max(nx,ny)
xlim=c(-nx/3,10)
ylim=c(0,top)
top <- max(nx,ny)
x <- list()
y <- list()
x.scale <- top/(nx+1)
y.scale <- top/(ny+1)
plot(NA,xlim=xlim,ylim=ylim,main=main,axes=FALSE,xlab="",ylab="")
for(i in 1:nx) {x[[i]] <- dia.rect(3,top-i*x.scale,x.names[i],cex=cex,...) }
 
for (j in 1:ny) {y[[j]] <- dia.rect(7,top-j*y.scale,y.names[j],cex=cex,...) }
for(i in 1:nx) {
  for (j in 1:ny) {
   dia.arrow(x[[i]]$right,y[[j]]$left,labels = beta[i,j],adj=4-j,cex=l.cex,...)
   }
} 
if(nx >1) {
  for (i in 2:nx) {
  for (k in 1:(i-1)) {dia.curved.arrow(x[[i]]$left,x[[k]]$left,x.matrix[i,k],scale=-(abs(i-k)),both=TRUE,dir="u",cex = l.cex,...)}
                      #dia.curve(x[[i]]$left,x[[k]]$left,x.matrix[i,k],scale=-(abs(i-k)))  } 
  } }
  
  if(ny>1) {for (i in 2:ny) {
  for (k in 1:(i-1)) {dia.curved.arrow(y[[i]]$right,y[[k]]$right,y.resid[i,k],scale=(abs(i-k)),dir="u",cex=l.cex, ...)} 
  }}
  for(i in 1:ny) {dia.self(y[[i]],side=3,scale=.2,... )}
 if(show) {text((10-nx/3)/2,0,paste("Unweighted matrix correlation  = ",round(sc$Ruw,digits)))}
}
```
Now try it on mod.2 (the traditional `lmCor` output) and then on mod.lm (the `lm` output)

```{r}
lmDiagram(mod.2, main="Regression from lmCor")
lmDiagram(mod.lm, main="Regression from lm")
```


# Code to speed up the graphics   

The original `setCor.diagram` function had a few tricks to make the graphics faster.  Rather than draw each shape and line separately we can not show show them until the last.

We pass parameters to the  various diagram functions saying not to draw, but rather collect the graphics to the end.  Although not a problem for small problems, it gets tedious for larger problems.

```{r}
lmDiagram <- setCor.diagram <- function(sc,main="Regression model",digits=2,show=FALSE,cex=1,l.cex=1,...) { 
if(missing(l.cex)) l.cex <- cex  

beta <- round(sc$coefficients,digits)
#two cases
#the normal case, the model comes from lmCor
#the 2nd case, the model comes from lm
if(length(class(sc)) ==1) {lm<- TRUE} else { lm <- FALSE}
if(!lm) {#case 1
if(rownames(beta)[1] %in% c("(Intercept)", "intercept*")) {intercept <- TRUE} else {intercept <- FALSE}
   
x.matrix <- round(sc$x.matrix,digits)

y.matrix <- round(sc$y.matrix,digits)
y.resid <- round(sc$resid,digits)
x.names <- rownames(beta)

if(intercept){ x.matrix <- x.matrix[-intercept,-intercept]
 x.names <- x.names[-intercept] 
 }
beta <- beta[-intercept,,drop=FALSE]
y.names <- colnames(beta)
} else { x.matrix <- round(cov(sc$model[-1]), digits=digits)  #case 2  from lmkl
       x.names <- colnames(sc$model)[-1]
       beta <- as.matrix(beta[-1]) #make the vector into a matrix
      y.names <-  colnames(sc$model)[1]
      }
      
nx <- length(x.names)
ny <- length(y.names) 
top <- max(nx,ny)
xlim=c(-nx/3,10)
ylim=c(0,top)
top <- max(nx,ny)
x <- list()
y <- list()
var.list <- arrow.list <-  curve.list <- self.list<- list()
x.scale <- top/(nx+1)
y.scale <- top/(ny+1)
plot(NA,xlim=xlim,ylim=ylim,main=main,axes=FALSE,xlab="",ylab="")
for(i in 1:nx) {x[[i]] <- dia.rect(3,top-i*x.scale,x.names[i],cex=cex,draw=TRUE,...) 
               var.list <- c(var.list,x.names[i],x[[i]])}
 
for (j in 1:ny) {y[[j]] <- dia.rect(7,top-j*y.scale,y.names[j],cex=cex,draw=TRUE,...)
                  var.list <- c(var.list,y.names[j],y[[j]]) }
for(i in 1:nx) {
  for (j in 1:ny) {
  
  d.arrow <-  dia.arrow(x[[i]]$right,y[[j]]$left,labels = beta[i,j],adj=4-j,cex=l.cex,draw=FALSE,...)
   arrow.list <- c(arrow.list,d.arrow)
   }
} 
if(nx >1) {
  for (i in 2:nx) {
  for (k in 1:(i-1)) {dca <- dia.curved.arrow(x[[i]]$left,x[[k]]$left,x.matrix[i,k],scale=-(abs(i-k)),both=TRUE,dir="u",cex = l.cex,draw=FALSE,...)
                     curve.list <- c(curve.list,dca)} #dia.curve(x[[i]]$left,x[[k]]$left,x.matrix[i,k],scale=-(abs(i-k)))  } 
  } }
  
  if(ny >1) {for (i in 2:ny) {
  for (k in 1:(i-1)) {dca <- dia.curved.arrow(y[[i]]$right,y[[k]]$right,y.resid[i,k],scale=(abs(i-k)),dir="u",cex=l.cex,draw=FALSE, ...)
     curve.list <- c(curve.list,dca)} 
  }}
  for(i in 1:ny) {d.self <- dia.self(y[[i]],side=3,scale=.2,draw=FALSE,... )
          self.list <- c(self.list,d.self)}
 if(show) {text((10-nx/3)/2,0,paste("Unweighted matrix correlation  = ",round(sc$Ruw,digits)))}
 
 
 #now draw all the arrows at once
 multi.rect(var.list,...)
 multi.arrow(arrow.list,...)
 multi.curved.arrow(curve.list,...)
 multi.self(self.list,...)
}
```


### test this one as well

```{r}
lmDiagram(mod.2, main="Regression from lmCor")
lmDiagram(mod.lm, main="Regression from lm")
diagram(mod.1)  #the call to the diagram function knows where to go
```

# Functions can redirect to other functions depending on the call

This concept is used in the `diagram` function which depending upon the object passed, will call one of some completely different functions.


```{r}
"diagram" <- function(fit,...) {
# figure out which psych function was called and then call the appropriate diagram
value <- NULL
omega <- bassAck <- extend <- extension <-  lavaan <- fa.reg <-  NULL    #weird requirement to avoid being global
#if(length(class(fit)) == 1) {if (inherits(fit,"lavaan")) fn <- "lavaan" } 


 #if(length(class(fit)) > 1)  {
    names <- cs(lavaan,fa,principal,iclust,omega,lavaan,bassAck, extension, extend, fa.reg, esem, setCor)
    value <- inherits(fit,names,which=TRUE)  # value <- class(x)[2]
    if(any(value > 1) ) { value <- names[which(value > 0)]} else {value <- "other"}
    
#     } else {value <- "other"}

switch(value,  #replaced a series of if statements to switch  12/23/18
fa  = {fa.diagram(fit,...) },
principal = {fa.diagram(fit,...) },
iclust = {iclust.diagram(fit,...)},
omega = {omega.diagram(fit,...)},
lavaan =  {lavaan.diagram(fit,...)}, 
bassAck = {bassAckward.diagram(fit,...)} ,
extend = {extension.diagram(fit ,...)},
extension ={extension.diagram(fit,...)},
fa.reg ={extension.diagram(fit,...)},
esem = {esem.diagram(fit,...)},
setCor = {setCor.diagram(fit,...)}
)
}

```

Lets try a few:

```{r}
f5 <- fa(bfi[1:25],5)
ic <- iclust(bfi[1:25])
om <- omega(ability,4, plot=FALSE)
#now draw these 3 output
diagram(f5)
diagram(ic)
diagram(om)
```