[R] pdIdent in smoothing regression model
Hi there,
I am reading the 2004 paper "Smoothing with mixed model software" in
Journal of Statistical Software, by Ngo and Wand. I tried to run
their first example in Section 2.1 using R but I had some problems.
Here is the code:
library(nlme)
fossil <- read.table("fossil.dat",header=T)
x <- fossil$age
y <- 10*fossil$strontium.ratio
knots <- seq(94,121,length=25)
n <- length(x)
X <- cbind(rep(1,n),x)
Z <- outer(x,knots,"-")
Z <- Z*(Z>0)
fit <- lme(y~-1+X,random=pdIdent(~-1+Z))
When I ran the code
fit <- lme(y~-1+X,random=pdIdent(~-1+Z))
I got an error message:
Error in getGroups.data.frame(dataMix, groups) : Invalid formula for groups
I was really puzzled. I asked Dr. Ngo and he said they did it in
S-plus but not R. Does anyone knows how to do it in R? Thanks!
Lei Liu
Associate Professor
Division of Biostatistics
Department of Public Health Sciences
University of Virginia School of Medicine
http://people.virginia.edu/~ll9f/
__
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.
[R] 2 questions on matrix manipulation in R
Hi there, I have two questions on matrix manipulation. For the first one, I want to calculate the product of each column of a matrix (say A) with another vector (say b). That is, if A has 5 columns (a1, a2, a3, a4, a5), I want to obtain a matrix with columns (a1*b, a2*b, aA3*b, a4*b, a5*b). Can I do it directly, without using "for" loop? For the second one, I have a matrix A of dimension 2 by n (say columns of a1, a2, ..., an), another matrix B of dimension 2 by 2. I want to obtain the vector with elements (t(a1) %*% B %*% a1, t(a2) %*% B %*% a2, ..., t(an) %*% B %*% an). Can I do it without using "for" loop? Thanks a lot! Lei Liu Associate Professor Division of Biostatistics Department of Public Health Sciences University of Virginia School of Medicine http://people.virginia.edu/~ll9f/ __ 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.
[R] multi-dimensional Gaussian quadrature
Hi there, Does anyone know if there is a package in R for multi-demensional Gaussian quadrature? I checked out the package "gaussquad" but it can only do 1D. Thanks! Lei Liu Associate Professor Division of Biostatistics Department of Public Health Sciences University of Virginia School of Medicine http://people.virginia.edu/~ll9f/ __ 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.
[R] How to install "adapt" package
Hi there, I want to install the "adapt" package, which is available at http://cran.r-project.org/src/contrib/Archive/adapt/. This package cannot be directly installed by "install packages" menu in R. So I downloaded the zip file into the hard drive. But I couldn't install it using "install packages from local zip files". I also tried to unzip these files in R library directory. But still I couldn't load the package in R. Is there any way to install this package? Thanks! Lei Liu Associate Professor Division of Biostatistics Department of Public Health Sciences University of Virginia School of Medicine http://people.virginia.edu/~ll9f/ __ 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.
[R] help on pspline in coxph
Hi there,
I have a question on how to extract the linear term in the penalized
spline in coxph. Here is a sample code:
n=100
set.seed(1)
x=runif(100)
f1 = cos(2*pi*x)
hazard = exp(f1)
T = 0
for (i in 1:100) {
T[i] = rexp(1,hazard[i])
}
C = runif(n)*4
cen = T<=C
y = T*(cen) + C*(1-cen)
data.tr=cbind(y,cen,x)
fit=coxph(Surv(data.tr[,1], data.tr[,2])~pspline(data.tr[,3]) )
If I use summary(fit), it will show the following results:
summary(fit)
Call:
coxph(formula = Surv(data.tr[, 1], data.tr[, 2]) ~ pspline(data.tr[,
3]))
n= 100
coef se(coef) se2 Chisq DF p
pspline(data.tr[, 3]), li 0.495 0.4370.437 1.28 1.00 2.6e-01
pspline(data.tr[, 3]), no 43.79 3.08 1.9e-09
exp(coef) exp(-coef) lower .95 upper .95
ps(data.tr[, 3])2 0.4404 2.270 0.08164 2.376
ps(data.tr[, 3])3 0.2065 4.842 0.01701 2.507
ps(data.tr[, 3])4 0.0951 10.512 0.00695 1.302
ps(data.tr[, 3])5 0.0493 20.274 0.00387 0.628
ps(data.tr[, 3])6 0.0280 35.741 0.00230 0.340
ps(data.tr[, 3])7 0.0192 52.068 0.00156 0.237
ps(data.tr[, 3])8 0.0219 45.605 0.00178 0.271
ps(data.tr[, 3])9 0.0473 21.156 0.00399 0.561
ps(data.tr[, 3])100.1432 6.983 0.01250 1.640
ps(data.tr[, 3])110.3936 2.541 0.03436 4.509
ps(data.tr[, 3])120.9449 1.058 0.0688512.969
ps(data.tr[, 3])132.2406 0.446 0.0764365.683
Iterations: 3 outer, 9 Newton-Raphson
Theta= 0.697
Degrees of freedom for terms= 4.1
Rsquare= 0.385 (max possible= 0.994 )
Likelihood ratio test= 48.6 on 4.08 df, p=7.74e-10
Wald test= 45.1 on 4.08 df, p=4.29e-09
My question is how to extract the linear coefficient (0.495) in
pspline(data.tr[, 3]). I tried coef(fit) but it fails to display this
term. Your help is greatly appreciated!
Lei Liu
Associate Professor
Division of Biostatistics and Epidemiology
Department of Public Health Sciences
University of Virginia School of Medicine
http://people.virginia.edu/~ll9f/
__
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.
[R] trajectory plot (growth curve)
Hi there, I want to make trajectory plots for data as follows: ID timey 1 1 1.4 1 2 2.0 1 3 2.5 2 1.5 2.3 2 4 4.5 2 5.5 1.6 2 6 2.0 ... That is, I will plot a growth curve for each subject ID, with y in the y axis, and time in the x axis. I would like to have all growth curves in the same plot. Is there a simple way in R to do it? Thanks a lot! Lei Liu Associate Professor Division of Biostatistics and Epidemiology Department of Public Health Sciences University of Virginia School of Medicine http://people.virginia.edu/~ll9f/ __ 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.
[R] a question in coxph
Hi there, I tried to fit a penalized spline for a continuous risk factor in recurrent events data. I found that I can include both pspline and frailty terms in coxph. So I use code like fit1 <- coxph(Surv(start, end, event)~pspline(age, df=0) + male + white +frailty(id), data.all) It yields results for both age spline and frailty variance. Next I tried to draw the plot for age effect. I used termplot function termplot(fit1, terms=c(1), se=TRUE) However, it gives me some error message: Error in dimnames(pred) <- list(dimnames(x)[[1]], termname) : length of 'dimnames' [2] not equal to array extent I don't know why. However, if I delete the frailty option, and use only fit2 <- coxph(Surv(start, end, event)~pspline(age, df=0) + male + white, data.all) I can generate the plot by termplot(fit2, terms=c(1), se=TRUE) So it seems to me that adding the frailty() option makes the "termplot" function fail to work? Do you know how to draw the plot for my function fit1? Thanks a lot! Lei At 09:53 AM 7/20/2009, you wrote: I would be willing to chair the session. Terry Therneau __ 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.
[R] package to do inverse probability weighting in longitudinal data
Hi there, I have a dataset from a longitudinal study with a lot of drop-out. I want to implement the inverse probability weighting method by Robins 1995 JASA paper "Analysis of semiparametric regression models for repeated outcomes in the presence of missing data". Does anyone know if there is a package to do it in R (or other software)? Thanks a lot! Lei __ 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.
[R] a question on matrix manipulation
Hi there, I have a question on manipulating a matrix. Say I have a matrix A with 3 rows. I want to generate a new matrix B with 3 duplicates of the first row of A, 2 duplicates of the second row, and 4 duplicates of the third row. So B is a matrix with 9 rows. Or more general, I want to generate (3, 2, 4, 5) duplicates of rows 1-4 of a matrix with 4 rows. Is there a simple function to do it? Thanks a lot! Lei __ 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.
