Nothing
# Find relationship between any two variables.
###############################
### The following function fits paired relations
###############################
## Input Arguments
# iVar: The first parameter.
# iResp: The second parameter.
# criterion: What model fit parameter is used for ranking? (Default: adj.R2)
# str: TRUE/FALSE Is this a 'strength' type model fitting? (whether to use ^2 and sqrt functions)
# modelNames: A vector of model names chosen from "SL", "Quad", "SQuad", "Exp", "Log", "nls","CP"
# cRes.all: A dataframe of all results, passed in by the netSEMm function
# cRes.best A dataframe of best results, passed in by the netSEMm function
# cRes.print A dataframe of results to print, passed in by the netSEMm function
# x A dataframe of data values, passed in by the netSEMm function
# nlsInits A vector of nls initialization coefficients, passed in by the netSEMm function
# nRes number of cells in the print variable, value passed in by the netSEm function ## (change netSEMm to netSEMp1 / see if it actually uses netSEMm)
## Return a list of the following items:
# "cRes.all": A dataframe of all results.
# "cRes.best": A dataframe of best results.
# "cRes.print": A dataframe of results to print
findrelation <- function(iVar,
iResp,
criterion = "adj.R2",
modelNames,
str = FALSE,
cRes.all,
cRes.best,
cRes.print,
x,
nlsInits,
nRes) {
#cat("\nRunning findrelation...\n")
#cat("\niVar is:",iVar)
#cat("\niResp is:",iResp,"\n")
## Save all adjusted R2 for all 8 functional forms
# Create zero vector with same dimension of modelNames
adj.R2 <- rep(0, length(modelNames))
## Save "lm" object for all 8 functional forms
model.res <- vector("list", length(modelNames))
Resp.v <- x[[iResp]] # Data of main endogenous variable
Var.v <- x[[iVar]] # Data of predictor
Resp <- colnames(x)[iResp] # Name of endogenous variable
Var <- colnames(x)[iVar] # Name of predictor
#cat("Regress Dep ~ Indep ", Resp, "~", Var, "\n")
lm4.flag <- TRUE
lm5.flag <- TRUE
lm6.flag <- TRUE
lm7.flag <- TRUE
lm9.flag <- TRUE
##--------------------------------------------------
## Functional Form 1: Linear
##--------------------------------------------------
Rel1 <- paste0(Resp, "~", Var)
lm1 <- do.call("lm", list(Rel1, data = as.name("x")))
adj.R2[1] <- summary(lm1)$adj.r.squared
model.res[[1]] <- lm1
cRes.all[[iVar, iResp, "SL"]] <- lm1
##--------------------------------------------------
## Functional Form 2: Quadratic
##--------------------------------------------------
Rel2 <- paste0(Resp, "~", Var, "+", "I(", Var, "^2)")
lm2 <- do.call("lm", list(Rel2, data = as.name("x")))
adj.R2[2] <- summary(lm2)$adj.r.squared
model.res[[2]] <- lm2
cRes.all[[iVar, iResp, "Quad"]] <- lm2
##--------------------------------------------------
## Functional Form 3: Quadratic (no linear term)
##--------------------------------------------------
Rel3 <- paste0(Resp, "~", "I(", Var, "^2)")
lm3 <- do.call("lm", list(Rel3, data = as.name("x")))
adj.R2[3] <- summary(lm3)$adj.r.squared
model.res[[3]] <- lm3
cRes.all[[iVar, iResp, "SQuad"]] <- lm3
##--------------------------------------------------
## Functional Form 4: Exponential
##--------------------------------------------------
# Test if lm4.flag is still true. To make sure exponential values are not infinity
lm4.flag <- all(exp(Var.v[!is.na(Var.v)]) != Inf)
if (lm4.flag) {
Rel4 <- paste0(Resp, "~", "exp(", Var, ")")
lm4 <- do.call("lm", list(Rel4, data = as.name("x")))
adj.R2[4] <- summary(lm4)$adj.r.squared
model.res[[4]] <- lm4
cRes.all[[iVar, iResp, "Exp"]] <- lm4
} else {
adj.R2[4] <- -Inf
model.res[[4]] <- NA
cRes.all[[iVar, iResp, "Exp"]] <- NA
}
##---------------------------------------------------
## Functional Form 5: Log
##---------------------------------------------------
lm5.flag <- all(Var.v > 0 & Var.v < Inf, na.rm = TRUE)
if (lm5.flag) {
Rel5 <- paste0(Resp, "~", "log(", Var, ")")
lm5 <- do.call("lm", list(Rel5, data = as.name("x")))
adj.R2[5] <- summary(lm5)$adj.r.squared
model.res[[5]] <- lm5
cRes.all[[iVar, iResp, "Log"]] <- lm5
} else {
adj.R2[5] <- -Inf
model.res[[5]] <- NA
cRes.all[[iVar, iResp, "Log"]] <- NA
}
# This is to be able to turn off these two functional forms
if (str) {
# Beginning of if(str)
##----------------------------------------------------
## Functional Form 6: SquareRoot (no linear term)
##----------------------------------------------------
lm6.flag <- all(Var.v > 0 & Var.v < Inf, na.rm = TRUE)
if (lm6.flag) {
Rel6 <- paste0(Resp, "~", "I(", Var, "^0.5)")
lm6 <- do.call("lm", list(Rel6, data = as.name("x")))
adj.R2[6] <- summary(lm6)$adj.r.squared
model.res[[6]] <- lm6
cRes.all[[iVar, iResp, "SQRoot"]] <- lm6
} else {
adj.R2[6] <- -Inf
model.res[[6]] <- NA
cRes.all[[iVar, iResp, "SQRoot"]] <- NA
}
##-------------------------------------------------------
## Functional Form 7: Inverse SquareRoot (no linear term)
##-------------------------------------------------------
lm7.flag <- all(Var.v > 0 & Var.v < Inf, na.rm = TRUE)
if (lm7.flag) {
Rel7 <- paste0(Resp, "~", "I(1/(", Var, "^0.5))")
lm7 <- do.call("lm", list(Rel7, data = as.name("x")))
adj.R2[7] <- summary(lm7)$adj.r.squared
model.res[[7]] <- lm7
cRes.all[[iVar, iResp, "ISQRoot"]] <- lm7
} else {
adj.R2[7] <- -Inf
model.res[[7]] <- NA
cRes.all[[iVar, iResp, "ISQRoot"]] <- NA
}
} # End of if(str)
##--------------------------------------------------------
## Functional Form 8: nls
##--------------------------------------------------------
## if(Resp == "G" & Var == "time") browser()
Rel8 <- paste0(Resp, " ~ ", "a1 + a2 * exp(a3 * ", Var, ")")
model8 <- NA
SSE <- NA
## Find the best (if any) nls models from all initial values
for (i in 1:nrow(nlsInits)) {
nls1 <-
tryCatch({
# Error Handling. Stores error message & output value.
do.call("nls", list(
Rel8,
data = as.name("x"),
start = list(
a1 = nlsInits[i, 1],
a2 = nlsInits[i, 2],
a3 = nlsInits[i, 3]
)
))
}, error = function(e)
e)
if (inherits(nls1, "nls")) {
# Check for "class" consistency
## If model8 is still NA, initialize it
if (!inherits(model8, "nls")) {
model8 <- nls1
SSE <- sum(residuals(nls1) ^ 2)
init <- nlsInits[i, ]
} else {
newSSE <- sum(residuals(nls1) ^ 2)
if (newSSE < SSE) {
model8 <- nls1
SSE <- newSSE
init <- nlsInits[i, ]
} else {
}
}
}
}
# This is needed to be able to drop out the 'str' functional forms
# and still have nls save into the cRes.all object correctly.
# It will be '8' if the 'str' offs occupy '6' and '7',
# but otherwise it will be '6'.
if (str) {
nlsNum <- 8
} else {
nlsNum <- 6
}
## Assign result
if (!inherits(model8, "nls")) {
adj.R2[nlsNum] <- -Inf
model.res[[nlsNum]] <- NA
cRes.all[[iVar, iResp, "nls"]] <- NA
} else {
R2s <-
nlsR2(model8, y = x[[Resp]], p = 3) ## Apply 'nlsR2' function
adj.R2[nlsNum] <- R2s$adjR2
model.res[[nlsNum]] <- model8
cRes.all[[iVar, iResp, "nls"]] <- model8
}
##--------------------------------------------------
## Functional Form 9: Change Point
##--------------------------------------------------
cpNum <- nlsNum + 1
lm9 <-
segmented.lm(
lm(Resp.v ~ Var.v, data = x),
seg.Z = ~ Var.v,
psi = NA,
control = seg.control(
K = 1,
fix.npsi = FALSE,
n.boot = 0,
it.max = 20
)
)
lm9.flag <-
all(summary(lm9)$adj.r.squared == summary(lm1)$adj.r.squared) # Evaluate logical condition
if (lm9.flag)
{
adj.R2[cpNum] <- -Inf
model.res[[cpNum]] <- NA
cRes.all[[iVar, iResp, "CP"]] <- NA
}
else
{
adj.R2[cpNum] <- summary(lm9)$adj.r.squared
model.res[[cpNum]] <- lm9
cRes.all[[iVar, iResp, "CP"]] <- lm9
}
## Choose the best model, if no model works, write -1.
## Populate the 'table' item in function return.
## Populate the 'bestModels' item in function return.
if (max(adj.R2, na.rm = TRUE) >= 0.001) {
step <- attributes(cRes.best)$Step[iResp] + 1
attributes(cRes.best)$Step[iVar] <-
min(step, attributes(cRes.best)$Step[iVar])
attributes(cRes.best)$diag.Step[iVar, iResp] <- step
nbest <- which.max(adj.R2) # Location of best model
cRes.best[iVar, iResp] <- modelNames[nbest]
## A row for the print table
table1r <- matrix(NA, nrow = 1, ncol = nRes)
colnames(table1r) <-
c(
"endogenous",
"Variable",
"Model",
"R-Sqr",
"adj-R-Sqr",
"Pval1",
"Pval2",
"Pval3"
)
table1r <- as.data.frame(table1r)
table1r[1, c("endogenous", "Variable")] = c(Resp, Var)
table1r[1, c("Model")] <- modelNames[nbest]
## if(modelNames[nbest] == "nls") browser()
## Best model from Resp ~ Var
finalModel <- model.res[[nbest]]
if (inherits(finalModel, "lm"))
{
model <- summary(finalModel)
table1r[1, c("R-Sqr", "adj-R-Sqr", "Pval1", "Pval2")] <-
c(model$r.sq, model$adj.r.sq, model$coef[1, 4], model$coef[2, 4])
if (modelNames[nbest] == "Quad" | modelNames[nbest] == "CP")
{
table1r[1, "Pval3"] <- model$coeff[3, 4]
}
else
{
table1r[1, "Pval3"] <- NA
}
}
else if (inherits(finalModel, "nls"))
{
table1r[1, c("R-Sqr", "adj-R-Sqr", "Pval1", "Pval2", "Pval3")] <-
c(R2s$R2, R2s$adjR2, NA, NA, NA)
}
cRes.print <- rbind(cRes.print, table1r)
} # Check if max adj.R2 >=0.001
else
# If not assign -1 for indicating no model
{
cRes.best[iVar, iResp] <- -1
}
return(list(cRes.all, cRes.best, cRes.print))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.