# demo/GetDerivs.R In dynr: Dynamic Models with Regime-Switching

```#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
# Illustrates how to get smoothed derivative estimates based on B-splines
# via wrapper function call to the fda package
# Author: Sy-Miin Chow
# The simulation model features:
# dx1(t)/dt = x2(t)
# dx2(t)/dt = eta*x1(t) + zeta*x2(t)
#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
require('dynr')

#Variables in this data set are:
#ID = ID of participants
#x = true scores
#theTimes = time indices

data("LinearOsc")

#Structure the one indicator for n individuals into a matrix with n columns
#If different individuals have different number of rows, consider doing
#this step separately for each individual

n = length(unique(LinearOsc\$ID)) #Number of subjects is 10
T = max(table(LinearOsc\$ID)) #Number of time points is 100
out2 = matrix(LinearOsc\$x,ncol=n,byrow=FALSE)
theTimes = LinearOsc\$theTimes[1:T]
norder = 6 #Order of Bsplines - usually 2 higher than roughPenaltyMax
roughPenaltyMax = 4 #penalization order
#  #lambdaLow, lambdaHi, lambdaBy = specify an interval of lambda (a positive smoothing parameter,
#                                 larger means more smoothing) to be tested from lambdaLow to
#                                 lambdaHi, as separated by lambdaBy
#isPlot = a binary flag on whether to plot the gcv values (0 = no, 1 = yes)
lambdaLow = 1e-10
lambdaHi = 2
lambdaBy = .1
isPlot = 1 #Whether to plot the average GCV values as a function of lambda values
matt = plotGCV(theTimes,norder,roughPenaltyMax,out2,lambdaLow, lambdaHi,lambdaBy,isPlot)

#Extract the lambda value that gives the minimum GCV
sp = matt[matt[,"GCV"] == min(matt[,"GCV"]),"lambda"]

#Extract smoothed level, first and second derivative estimates at the lambda value selected above
x = getdx(theTimes,norder,roughPenaltyMax,sp,out2,0)[] #Smoothed level
dx = getdx(theTimes,norder,roughPenaltyMax,sp,out2,1)[] #Smoothed 1st derivs
d2x = getdx(theTimes,norder,roughPenaltyMax,sp,out2,2)[] #Smoothed 2nd derivs

#Put level and derivative estimates into a data frame
dxall = data.frame(time = rep(theTimes,n),
x = matrix(x,ncol=1,byrow=FALSE),
dx = matrix(dx,ncol=1,byrow=FALSE),
d2x = matrix(d2x,ncol=1,byrow=FALSE))

g = lm(d2x~x+dx-1,data=dxall)

#Component + residuals plot to show the association between smoothed d2x and smoothed x
#after partialling out the effect of smoothed dx
oldpar <- par(mgp=c(2.5,0.5,0))
car::crPlots(g,terms=~x,
main=c(""),layout=c(1,1),
cex.lab=1.3,cex.axis=1.2,
xlab=expression(hat(eta)[i](t)),
ylab=expression(paste("Component+Residuals ", "  ",d^2,hat(eta)[i](t)/dt^2))
)
par(oldpar)

#Component + residuals plot to show the association between smoothed d2x and smoothed dx
#after partialling out the effect of smoothed x
oldpar <- par(mgp=c(2.5,0.5,0))
car::crPlots(g,terms=~dx,
main=c(""),layout=c(1,1),
cex.lab=1.3,cex.axis=1.2,
xlab=expression(paste(d, hat(eta)[i](t)/dt)),
ylab=expression(paste("Component+Residuals ", "  ",d^2,hat(eta)[i](t)/dt^2))
)
par(oldpar)

# ---- Plot of simple slopes and region of significance ----
g2 = lm(d2x~x+dx+x:dx-1,data=dxall) #Adding an interaction term to illustrate some
# functions for probing interaction effects
summary(g2) #In this case the data were generated without any interaction effect.
#With larger T or n, sometimes spurious interaction effects may be detected
#due to shared variability between x and dx.
theta_plot(g2, predictor = "x", moderator = "dx",
alpha = .05, jn = T, title0=" ",
predictorLab = "x", moderatorLab = "dx")

# ---- Phase portrait ----
Osc <- function(t, y, parameters) {
dy <- numeric(2)
dy <- y
dy <- parameters*y+parameters*dy
return(list(dy))
}

param <- coef(g)
dynr.flowField(Osc, xlim = c(-3, 3),
ylim = c(-3, 3),
xlab="x", ylab="dx/dt",
main=paste0("Oscillator model"),
cex.main=2,
parameters = param,
points = 15, add = FALSE,
col="blue",
arrow.type="proportional",