3
$\begingroup$

I am trying to create a predictive model using a relationship found by the lm.circular function from the circular package. The function, however, only provides statistics of the models fit and the intercept/coefficients, and I have been struggling to put together the equation of the linear model using these outputs. How can I take the outputs from this function and develop an equation to predict the y variable?

I've done some reading inJammalamadaka, S. Rao and SenGupta, A. (2001), Downs and Mardia (2002) and have included an example below taken from Pewsey, Neuhauser and Graeme (2013), but am still unsure of how to construct those relationships from what R is giving me.

library(circular)

psideg <- c(356,97,211,232,343,292,157,302,335,302,324,85,324,340,157,238,254,146,232,122,329)

thetadeg <- c(119,162,221,259,270,29,97,292,40,313,94,45,47,108,221,270,119,248,270,45,23)

cpsirad <- circular(psideg*2*pi/360); cthetarad <- circular(thetadeg*2*pi/360)

circ.lm <- lm.circular(type = "c-c", y = cthetarad, x = cpsirad, order = 1)

plot.default(cthetarad, cpsirad, xlab = "cthetarad", ylab = "cpsirad", xlim = c(0, 7), ylim = c(-2, 7))

circ.lm$rho
#>          [,1]
#> [1,] 0.502289

circ.lm$coefficients
#>                    [,1]      [,2]
#> (Intercept)  0.01441183 0.1145101
#> cos.x        0.33465348 0.5494121
#> sin.x       -0.07478683 0.4821004
$\endgroup$
0
$\begingroup$

Using edits and additions from Ulric Lund, one of the circular package co-authors, the following script utilizes the model coefficients to recreate the model manually for the purposes of building a predictive relationship.

## Checking how to get fitted values
### Get model coefficients
( cf <- circ.lm$coefficients )
#>                    [,1]      [,2]
#> (Intercept)  0.01441183 0.1145101
#> cos.x        0.33465348 0.5494121
#> sin.x       -0.07478683 0.4821004

### Model calculations provided by package co-author Ulric Lund
costheta <- cf[1,1]+cf[2,1]*cos(psirad)+cf[3,1]*sin(psirad)
sintheta <- cf[1,2]+cf[2,2]*cos(psirad)+cf[3,2]*sin(psirad)
# Quadrant-correct version of atan:
calcfits <- atan2(sintheta, costheta)

The calcfits is the calculated values of the manual model and the diffits shows the fit between models. The difference of 6.28319 demonstrates the circular nature of the data and can be accounted for by adding 2pi because of the 1:1 relationship as visualized in the plot outputs.

### Difference between model and hand-calculated fitted values
( difffits <- round(circ.lm$fitted-calcfits,5) )
#> Circular Data: 
#> Type = angles 
#> Units = radians 
#> Template = none 
#> Modulo = asis 
#> Zero = 0 
#> Rotation = counter 
#>       1       2       3       4       5       6       7       8       9 
#> 0.00000 0.00000 6.28319 6.28319 0.00000 6.28319 6.28319 6.28319 0.00000 
#>      10      11      12      13      14      15      16      17      18 
#> 6.28319 0.00000 0.00000 0.00000 0.00000 6.28319 6.28319 6.28319 6.28319 
#>      19      20      21 
#> 6.28319 0.00000 0.00000

## Making a fitted line plot from the model results
plot.default(cpsirad, cthetarad, xlab = "cpsirad", ylab = "cthetarad", xlim = c(0, 6.5), ylim = c(0, 6.5))
lines.default(cpsirad[order(cpsirad)], circ.lm$fitted[order(cpsirad)])

enter image description here

## ... same data but adjusted for 2*pi
thetaPlot  <- ifelse(cpsirad < 5.5, cthetarad, cthetarad + 2*pi)
fittedPlot <- ifelse(cpsirad < 5.5, circ.lm$fitted, circ.lm$fitted + 2*pi)
plot.default(cpsirad, thetaPlot, xlab = "cpsirad", ylab = "cthetarad")
lines.default(cpsirad[order(cpsirad)], fittedPlot[order(cpsirad)])

enter image description here

Other notes include a error in the original plot code that had theta and psi switched, and further information on using p values. The output p values from the package test whether a higher order trigonometric polynomial is warranted. The null hypotheses in this case is that the current model (with only first order terms) is sufficient. The p-values are not small, so we don't reject these null hypotheses, and so we treat the estimated model as the final model.

circ.lm$p.values
#>            [,1]      [,2]
#> [1,] 0.09971087 0.9474847

#Rho output can be used similarly to r^2
circ.lm$rho
#>          [,1]
#> [1,] 0.502289
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.