Dylan Beaudette wrote:
> On Mon, Dec 1, 2008 at 5:48 AM, Frank E Harrell Jr
> <[EMAIL PROTECTED]> wrote:
>> David Winsemius wrote:
>>> On Nov 30, 2008, at 10:23 PM, Dylan Beaudette wrote:
>>>
>>>> Hi, I am using the rcs() function in the Design library to model
>>>> non-linearity that is not well characterized by an otherwise
>>>> mechanistic function. I am able to make the model 'available' to
>>>> others through the excellent nomogram() function and the set of tables
>>>> that it can create. However, I would like to present the model in an
>>>> 'expanded' format-- probably what rcspline.restate() or latex.Design()
>>>> produce on a model fit object.
>>>>
>>>> Here is how the model was fit:
>>>>
>>>> fit.ols <- ols( log(k) ~ (rcs(activity) * (log(conc) + sar)) +
>>>> (rcs(sand) * (log(conc) + sar)), data=sm.clean, x=TRUE, y=TRUE)
>>>>
>>>> Here is how I am accessing the 'expanded' format of the model structure:
>>>>
>>>> options(digits=3)
>>>> latex(fit.ols, file='fit_rcs.tex')
>>>>
>>>> The output contains several notation elements that I am not familiar
>>>> with:
>>>>
>>>> 1. x_{+} --> it seems that this represent a term that should be set
>>>> to 0, when x is 0?
>>> It is set to zero when the term inside the cubic is less than zero. See
>>> pages 20-21 of Harrell's book where the basis functions are described
>>> and illustrated.
>>>> i.e. the entire expression −453(activity − 0.842)_{+}^{3} = 0 when
>>>> 'activity' = 0 ??
>>> .... whenever (activity − 0.842) < 0
>>>
>>>
>>>
>>>> 2. the '!x' found in :
>>>>
>>>> +log(conc) [ −0.0118sand + 9.58
>>>> ! ×
>>>> !10−6 (sand − 11.6)
>>> My guess is that this is 9.58 x 10^-6
>>>> − 0.000128(sand − 37.5)
>>>> +0.00045(sand − 47.2)
>>>> − 0.000350(sand − 51)
>>>> + 1.86
>>>> ! ×
>>>> !10−5 (sand − 69.8) ]
>>> I don't see anything like that in Harrell's text and I am wondering if a
>>> different character is getting rendering incorrectly. The only time you
>>> see it is when the exponent is below -4.
>>>
>>>
>>>>
>>>>
>>>> .... what exactly does that mean?
>>>>
>>>>
>>>> An image version of the equation in question is attached.
>>>>
>>>> Any input would be greatly appreciated!
>>>>
>>>> Cheers,
>>>>
>>>> Dylan
>> Thanks for the good reply David.
>>
>> There is a bug in one of the Hmisc functions regarding the ! spacing
>> command. If you get any LaTeX syntax errors I can send you
>> the correction.
>>
>> Frank
>
> I didn't get any show-stopping errors when compiling the latex code
> used to generate the previously attached equation. Looking through the
> LATEX log file I don't see anything related to the equation and odd
> '!x' stuff.
>
> HOWEVER: I think that this is where the error is occuring:
>
> \\!\times\\!10^{-6}
>
> * should be *
>
> \! \times \! 10^{-6}
>
> What do you think Frank? Altering the above syntax results in a
> normal-looking presentation.
Yes you can remove the extra \ or source( ) in the attached new version
of one of the functions that will be in an upcoming release of Hmisc.
Frank
>
> Cheers,
>
> Dylan
>
>
rcspline.restate <- function(knots, coef, type=c("ordinary","integral"),
x="X", lx=nchar(x),norm=2,
columns=65, before="& &", after="\\",
begin="", nbegin=0,
digits=max(8,.Options$digits))
{
type <- match.arg(type)
k <- length(knots)
if(k<3)
stop("must have >=3 knots in a restricted cubic spline")
p <- length(coef)
if(p == k)
{
Intc <- coef[1]
coef <- coef[-1]
p <- p-1
}
else Intc <- 0
if(k-1 != p)
stop("coef must be of length # knots - 1")
knotnk <- knots[k]
knotnk1 <- knots[k-1]
knot1 <- knots[1]
kd <- if(norm==0) 1 else if(norm==1)(knotnk-knotnk1)^3 else (knotnk-knot1)^2
coef[-1] <- coef[-1]/kd
d <- c(0, knots-knotnk)[1:p]
coefk <- sum(coef*d)/(knotnk-knotnk1)
d <- c(0, knots-knotnk1)[1:p]
coefk1 <- sum(coef*d)/(knotnk1-knotnk)
if(!length(names(coef)))
names(coef) <- paste(x,1:length(coef),sep="")
coef <- c(coef, coefk, coefk1)
names(coef)[k] <- "1st restricted coef"
names(coef)[k+1] <- "2nd restricted coef"
if(type=="integral")
coef <- c(.5*coef[1],.25*coef[-1])
cof <- format.sep(coef, digits)
kn <- format.sep(-knots, digits)
if(Intc!=0)
{
txt <- txt2 <- format.sep(Intc, digits)
if(type=="integral")
{
txt <- paste(txt, "* x")
txt2 <- paste(txt2, '*', x)
}
if(coef[1]>=0)
{
txt <- paste(txt, "+");
txt2 <- paste(txt2, '+')
}
}
else txt <- txt2 <- ""
if(cof[1]!=0)
{
txt <- paste(txt, cof[1],
if(type=="ordinary")"* x"
else "* x^2",
sep="")
txt2 <- paste(txt2, cof[1],
if(type=="ordinary") paste("*",x)
else paste("*",x,"^2"),
sep="")
}
for(i in 2:(p+2))
{
nam <- paste("pmax(x",
if(knots[i-1]<0) "+"
else NULL,
if(knots[i-1]!=0) kn[i-1]
else NULL,
",0)^",
if(type=="ordinary")"3"
else "4",
sep="")
nam2 <- paste("pmax(",x,
if(knots[i-1]<0) "+"
else NULL,
if(knots[i-1]!=0) kn[i-1]
else NULL,
",0)^",
if(type=="ordinary")"3"
else "4",
sep="")
z <- paste(if(coef[i]>0 & (i>2 | coef[1]!=0 | Intc!=0)) "+"
else NULL,
cof[i], "*", nam, sep="")
z2 <- paste(if(coef[i]>0 & (i>2 | coef[1]!=0 | Intc!=0)) "+"
else NULL,
cof[i], "*", nam2, sep="")
txt <- paste(txt , z, sep="")
txt2<- paste(txt2, z2, sep="")
}
func <- parse(text=paste('function(x)', txt))
cof <- format.sep(coef, digits)
kn <- format.sep(-knots, digits)
lcof <- nchar(cof)
cof <- latexSN(cof)
cur <- begin; colcnt <- nbegin; tex <- NULL
if(Intc!=0)
{
fint <- format.sep(Intc, digits)
if(type=="integral")
{
fint <- paste(fint, x)
colcnt <- colcnt+2
}
cur <- paste(cur, fint, sep="")
colcnt <- colcnt + nchar(fint)
if(coef[1]>0)
{
cur <- paste(cur, " + ", sep="");
colcnt <- colcnt+3
}
}
if(coef[1]!=0)
{
sp <- if(substring.location(cof[1],"times")$first > 0) "\\:"
else NULL
cur <- paste(cur, cof[1], sp, x,
if(type=="integral") "^2",
sep="")
##\:=medium space in LaTeX
colcnt <- colcnt+lcof[1]+lx+(type=="integral")
}
tex.names <- character(p+2)
size <- lx+lcof[-1]+nchar(kn)+3
for(i in 2:(p+2))
{
nam <- paste("(", x,
if(knots[i-1]<0) "+"
else NULL,
if(knots[i-1]!=0) kn[i-1]
else NULL,
")_{+}^{",
if(type=="ordinary")"3}"
else "4}",
sep="")
q <- paste(if(coef[i]>0 & (i>2 | coef[1]!=0 | Intc!=0)) "+"
else NULL,
cof[i], nam, sep="")
n <- size[i-1]
if(colcnt+n > columns)
{
tex <- c(tex, cur)
cur <- ""
colcnt <- 0
}
cur <- paste(cur, q, sep="")
colcnt <- colcnt+n
}
tex <- c(tex, cur)
tex <- paste(before, tex, after)
if(Intc!=0) coef <- c(Intercept=Intc, coef)
attr(coef, "knots") <- knots
attr(coef, "function") <- func
attr(coef, "function.text") <- txt2
attr(coef, "latex") <- tex
names(colcnt) <- NULL
attr(coef, "columns.used") <- colcnt
coef
}
rcsplineFunction <- function(knots, coef=numeric(0), norm=2)
{
k <- length(knots)
kd <- if(norm==0) 1 else if(norm==1) knots[k]-knots[k-1] else
(knots[k]-knots[1])^.66666666666666666666666
f <- function(x, knots, coef, kd)
{
k <- length(knots)
knotnk <- knots[k]
knotnk1 <- knots[k-1]
knot1 <- knots[1]
if(length(coef) < k) coef <- c(0, coef)
y <- coef[1] + coef[2]*x
for(j in 1:(k-2))
y <- y +
coef[j+2]*(pmax((x - knots[j])/kd, 0)^3 +
((knotnk1 - knots[j]) * pmax((x - knotnk)/kd, 0)^3 -
(knotnk - knots[j]) * (pmax((x - knotnk1)/kd, 0)^3))/
(knotnk - knotnk1))
y
}
formals(f) <- list(x=numeric(0), knots=knots, coef=coef, kd=kd)
f
}
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
