1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
|
##rcspline.eval - function to create design matrix for restricted cubic
## spline function of Stone & Koo, given an input vector and optionally
## a vector of knots. If knots are not given, knots are set using
## default algorithm. If the number of knots is not given, 5 are used.
## Terms are normalized by (outer-inner knot)^2.
## Can optionally return antiderivative of spline functions if
## type="integral".
## norm=0 : no normalization of constructed variables
## norm=1 : divide by cube of difference in last 2 knots
## makes all variables unitless
## norm=2 : (default) divide by square of difference in outer knots
## makes all variables in original units of x
##
## Returns:
## x - design matrix for derived spline variables
## (includes original x in first column if inclx=T or
## type="integral")
## attribute knots - input or derived vector of knots
## If knots.only=T, returns instead the vector of estimated or given
## knots.
## If rpm is not null, replaces missing x with rpm before evaluating
## but after estimating knots.
##
## F. Harrell 13 Feb 90
## Modified 28 Mar 90 - improved default knot computation
## 22 Aug 90 - put knots as attribute, return matrix
## 20 Sep 90 - added knots.only argument
## 16 Oct 90 - added rpm argument
## 11 Dec 91 - added type argument
## 27 Dec 91 - added norm argument
## 26 Jun 93 - added evasive action if <3 knots
rcspline.eval <- function(x,knots=NULL,nk=5,inclx=FALSE,knots.only=FALSE,
type="ordinary",norm=2, rpm=NULL, pc=FALSE)
{
if(!length(knots))
{
xx <- x[!is.na(x)]
n <- length(xx)
if(n<6)
stop('fewer than 6 non-missing observations with knots omitted')
if(nk<3)
stop('nk must be >= 3')
outer <- if(nk > 3) .05 else .1
if(nk>6) outer <- .025
knots <- quantile(xx,seq(outer,1.0-outer,length=nk))
if(length(unique(knots))<3)
{
knots <- quantile(xx,seq(outer,1.0-outer,length=2*nk))
if((nu <- length(unique(knots)))<3) {
cat("Fewer than 3 unique knots. Frequency table of variable:\n")
print(table(xx))
stop()
}
warning(paste("could not obtain",nk,"knots with default algorithm.\n",
"Used alternate algorithm to obtain",
nu,"knots"))
}
if(n<100)
{
xx <- sort(xx)
knots[1]<-xx[5]
knots[nk]<-xx[n-4]
}
}
knots <- sort(unique(knots))
nk <- length(knots)
if(nk<3)
{
cat("fewer than 3 unique knots. Frequency table of variable:\n")
print(table(x))
stop()
}
if(knots.only) return(knots)
if(length(rpm)) x[is.na(x)] <- rpm
xx <- matrix(1.1,length(x),nk-2)
knot1 <- knots[1]
knotnk <- knots[nk]
knotnk1 <- knots[nk-1]
kd <- if(norm==0) 1
else if(norm==1) knotnk-knotnk1
else (knotnk-knot1)^(2/3)
power <- if(type=="integral")4 else 3
for(j in 1:(nk-2))
{
xx[,j]<-pmax((x-knots[j])/kd,0)^power +
((knotnk1-knots[j])*pmax((x-knotnk)/kd,0)^power -
(knotnk-knots[j])*(pmax((x-knotnk1)/kd,0)^power))/
(knotnk-knotnk1)
}
if(power==4) xx <- cbind(x, x*x/2, xx*kd/4)
else if(inclx) xx <- cbind(x, xx)
if(!.R.) storage.mode(xx) <- 'single'
if(pc)
{
p <- prcomp(xx, scale=TRUE, center=TRUE)
pcparms <- p[c('center','scale','rotation')]
xx <- p$x
attr(xx, 'pcparms') <- pcparms
}
attr(xx, 'knots') <- knots
xx
}
|
