0

Following the very helpful discussion here:

Export fitted regression splines (constructed by 'bs' or 'ns') as piecewise polynomials

I would like to use the R package SplinesUtils to export estimates of the piecewise polynomials for regression models that include interaction terms with the spline.

Here is the model I am fitting:

scholl_P9h =lme(number~bs(radius, knots=c(60), degree=2)*genotype, data = P9_300,random=list(mouse= ~1, cell=~1),na.action = "na.omit")

scholl_P9h
Linear mixed-effects model fit by REML
  Data: P9_300 
  Log-restricted-likelihood: -2379.48
  Fixed: number ~ bs(radius, knots = c(60), degree = 2) * genotype 
                                      (Intercept)            bs(radius, knots = c(60), degree = 2)1 
                                        0.5193445                                        25.9796546 
           bs(radius, knots = c(60), degree = 2)2            bs(radius, knots = c(60), degree = 2)3 
                                        4.4489480                                        -3.1056313 
                                       genotypeWT bs(radius, knots = c(60), degree = 2)1:genotypeWT 
                                        1.5858730                                        -5.7978228 
bs(radius, knots = c(60), degree = 2)2:genotypeWT bs(radius, knots = c(60), degree = 2)3:genotypeWT 
                                       -4.3320939                                         0.2900458 

Random effects:
 Formula: ~1 | mouse
        (Intercept)
StdDev:    2.377891

 Formula: ~1 | cell %in% mouse
        (Intercept) Residual
StdDev:    5.190573 4.244785

Number of Observations: 812
Number of Groups: 
          mouse cell %in% mouse 
              8              37

I am able to get the coefficients for the degree 2 piecewise polynomial both before and after the knot for the group with genotype = 0

> ans1 <- RegSplineAsPiecePoly(scholl_P9h, "bs(radius, knots = c(60), degree = 2)")
> ans1
2 piecewise polynomials of degree 2 are constructed!
Use 'summary' to export all of them.
The first 2 are printed below.
4.1e-15 + 0.866 * (x + 0) - 0.00841 * (x + 0) ^ 2
21.7 - 0.144 * (x - 60) + 0.000168 * (x - 60) ^ 2

But I can't do the same for the group with genotype = 1.

> ans2 <- RegSplineAsPiecePoly(scholl_P9h, "bs(radius, knots = c(60), degree = 2)*genotype")
Error: 
  Required SplineTerm not found! Available terms are:
    * bs(radius, knots = c(60), degree = 2)
    * genotype
    * bs(radius, knots = c(60), degree = 2):genotype
> ans2 <- RegSplineAsPiecePoly(scholl_P9h, "bs(radius, knots = c(60), degree = 2):genotype")
Error in predvars[[2L + pos]] : subscript out of bounds

Any suggestions?

farzad
  • 1

0 Answers0