Re: [R] 3D plot and interactive PDFs

This does not answer exactly to your question, anyway..: If you are planning to use latex, the package movie15 allows to include media files in your document (to be processed via pdflatex) vito Bruno C. wrote: With version 8 of acrobat reader, it is now possible to have 3D in PDf

Re: [R] Latex \ell symbol in plotmath

Dear Mario, I don't know whether the $\ell$ symbol is available .. However you can use the LaTeX psfrag/pdfrag packages to convert tags in latex symbols.. hope this helps, vito Mario dos Reis wrote: Is it possible to use the '\ell' (i.e. the log likelihood) in plots? I've been browsing the

Re: [R] Testing additive nonparametric model

You can use the LRT (although I think that it assumes the df to be fixed). For instance the package mgcv by Simon Wood has an anova method to compare models fitted by the relevant gam() function, and the print.summary() itself returns such information.. best, vito set.seed(123) n-100 sig-2 x0

Re: [R] Any R function for self-controlled case series method /effect absorption?

Dear Jari The problem is to build the dataset to apply the conditional logit model. However, as far as I know, no R function exists. BTW if you are dealing with time series of pollution and health, the following two papers might be of interest of you: It appears that the time series approach

[R] difference in using with() and the data argument in glm call

Dear all, I am dealing with the following (apparently simple problem): For some reasons I am interested in passing variables from a dataframe to a specific environment, and in fitting a standard glm: dati-data.frame(y=rnorm(10),x1=runif(10),x2=runif(10)) KK-new.env() for(i in 1:ncol(dati))

Re: [R] Question about managing searching path

Dear Tong I think on.exit() makes the job..Namely: attach(Yourdata) on.exit(detach(YourData)) vito Tong Wang wrote: Hi all, I'm having sometrouble with managing the seach path, in a function , I need to attach some data set at the begining and detach them at the end, say,

Re: [R] Constrained OLS regression

In addition to Dimitris's approach, probably the following is more straightforward..(the idea is the same, but implementation is simpler; you do not need starting values, for instance..) Given the linear predictor lp: b0+b1X1+b2X2 as b2=1-b1 the lp becomes: b0+b1X1+(1-b1)X2 =

Re: [R] How to generating diagnal blocks ?

If I remember well, there should be a package including the function bdiag() making the job..but I am not able to remember its name.. However a quick search via RSiteSearch(bdiag) yields http://finzi.psych.upenn.edu/R/Rhelp02a/archive/40393.html Hope this helps you, vito Tong Wang wrote:

Re: [R] unexpected result in glm (family=poisson) for data with an only zero response in one factor

Dear Antonin, It is a statistical problem: the well-known monotone likelihood. In this case ML estimate does not exist (or equals infinity) and Wald approximations (ob which SE are based) does not hold. However LRT is valid and provides reliable results. As far as I know, the only software

Re: [R] clogit question

Dear Dan, I think you need more (theorical) background here.. clogit() in package survival performs conditional logistic regression where you have several groups (the strata, the matched sets). There is an intercept for each stratum in the model, but you do not obtain them since estimation is

Re: [R] matrix indexing

Dear tom, is the following what you are looking for? a=matrix(runif(9),3,3) a [,1] [,2] [,3] [1,] 0.9484247 0.9765431 0.6169739 [2,] 0.8423545 0.3137295 0.4031847 [3,] 0.6724235 0.1076373 0.2356923 b-matrix(sample(c(TRUE,FALSE),size=9,replace=TRUE),3,3) b [,1]

[R] a strange problem with integrate()

Dear all, I am stuck on the following problem with integrate(). I have been out of luck using RSiteSearch().. My function is g2-function(b,theta,xi,yi,sigma2){ xi-cbind(1,xi) eta-drop(xi%*%theta) num-exp((eta + rep(b,length(eta)))*yi) den- 1 + exp(eta +

Re: [R] how to use the basis matrix of ns in R? really confused by multi-dim spline filtering?

Dear Micheal, the output of the ns function in R is basis matrix, but then Yes you are right, the output of the ns(x, df) is the basis matrix of a natural cubic spline with df degrees of freedom. See ?ns (in package splines) on how to specify df or knots or .. Fitting y~ns(x,df) yields a

Re: [R] warnings in glm (logistic regression)

Hi, very likely your data exhibit quasi-separation which cause (log)Lik to be monotone and thus ML estimate do not exist. However you can rely on point estimate and use LRT to test for its significance. Or Better: have a look to brlr or logistf packages which bypass the monotone-likelihood

Re: [R] Breakpoints for multiple variables using Segmented

Dear Matthew, Currently segmented() performs multiple (say L1) estimation of breakpoints in GLM, namely: 1)L breakpoints for the same variable x, e.g.: segmented(obj.glm, Z=x, psi=c(psi1,psi2,psi3)) 2)L breakpoints for L explanatory `segmented' variables, e.g.: segmented(obj.glm,

Re: [R] AIC in lmer

Hi, my reply just concerns the usage of AIC in mixed models and not the lmer package. The standard AIC is actually unconditional. Vaida and Blanchard (2003, Proceeding 19 IWSM,101-105) discuss that a conditional version should be more appropriate in a mixed framework. I don't whether the

Re: [R] Robustness of Segmented Regression Contributed by Muggeo

Hi, sorry for my delay.. In addition to valuable Achim's comments. As Achim said, you can try different starting values to assess how the final solution depends on them. Then select one having the best logLik (or the minimum RSS). Everybody dealing with nonlinear models knows that the logLik

Re: [R] grubbs.test

Dear Dave, I do not know the grubbs.test (is it a function, where can I find it?) and probably n=6 data points are really few.. Having said that, what do you mean as outlier? If you mean deviation from the estimated mean (of previous data), you might have a look to the strucchange

Re: [R] Fitting mixed proportional odds model in R?

Zoran Loncarevic wrote: Is there a way to fit mixed proportional odds models in R? As far as I know, no. (anyway have a look to J Lindsey's packages, I don't know) However MIXOR and friends at http://tigger.uic.edu/~hedeker/mix.html (standalone programs running on Windows systems ) can fit mixed

Re: [R] Basic stratification calculations

Hi, if I understand correctly, tapply() is your friend here, vito [EMAIL PROTECTED] wrote: Hi. I'm a student at SFU in Canada. The basic thing I want to do is calculate means of different strata. I have 2 vectors. One has the values I want to take the means from, the other is the four strata I am

Re: [R] Negative binomial regression for count data

Hi, I do not know the article. Notice that an excess of zeroes can lead to (spurious) overdispersion in data, therefore you should decide whether assuming a zip ( zero excess coming from a mixture) or a negBin (zero execess due to overdispersion) model. Of course some likelihood based criteria

R: [R] partial linear model

-glm(y~X+_someKnownFunction(th)_+..) o$dev } #search the optimum ob-optimize(fn,.. th1-ob$minimum #(or ob$maximum) o-glm(y~X+_someKnownFunction(th1)_+..) #fit the model assuming th=th1 *known* Hope this helps, vito muggeo - Original Message - From: Jin Shusong [EMAIL PROTECTED] To: R

R: [R] time dependency of Cox regression

Hi, in order to fit Cox model with time-dependent coeff, you have to restruct your dataframe. For instance you can use the counting process formulation (start,stop,status). Some years ago I wrote a function (reCox() below) to make the job. It seems to work if there are not ties, but if there are

R: [R] loess significance

Dear Aurélie, I think that for *fixed* (i.e. assumed known) amount of smoothing, you can use a simple LRT by comparing the two candidate models. BTW, have a look to the mgcv or gam packages for a general model-based approach. - Original Message - From: Aurélie Coulon [EMAIL PROTECTED]

R: [R] Indexing lists

Hi, Probably it should be useful to obtain two results which you can combine according to what you need. Here a possible solution: a-b-data.frame(matrix(runif(30),10,3)) d-list(a,b) #your list #extract the 1st column. The output is a nrow(a)-by-length(d) matrix sapply(d,function(x)x[,1])

[R] an integration question

Dear all, I'm stuck on a problem concerning integration..Results from the analytical expression and numerical approximation (as returned by integrate()) do not match. It probably depends on some error of mine, so apologizes for this off-topic question. I'm interested in computing the integral of

R: [R] solving for a transition point in a piecewise nonlinear model

Dear Robert, In general it may be difficult to estimate a model with generic (possibly nonlinear) functions before/after the changepoint to be estimated too. However if you are willing to make some restrinctions on your F1(.) and F2(.), you could semplify the problem.. For instance, have a look

R: [R] slope estimations of teeth like data

Dear Petr, Probably I don't understand exactly what you are looking for. However your plot(x,c(y,z)) suggests a broken-line model for the response c(y,x) versus the variables x. Therefore you could estimate a segmented model to obtain (different) slope (and breakpoint) estimates. See the package

R: [R] column sorting a matrix with indices returned

A possible solution is using apply + order: apply(x,2,order) [,1] [,2] [,3] [1,]323 [2,]132 [3,]211 apply(x,2,function(x)x[order(x)]) [,1] [,2] [,3] [1,]235 [2,]347 [3,]468 best, vito - Original Message

R: [R] absolute value

Are you looking for Re() and friends? Toy examples: abs(Re(3+4i)) [1] 3 abs(Re(-3+4i)) [1] 3 abs(Re(3)) [1] 3 see ?complex for further details on complex numbers in R vito - Original Message - From: Fred J. [EMAIL PROTECTED] To: r help [EMAIL PROTECTED] Sent: Friday, April 30,

[R] random effects structure in lme

Dear all, Some days ago I posted the same message below, but unfortunately with no reply. So, apologizes for this my re-sending, but I hope there is now someone on-line that can help me In using the lme() function from the nlme package, I would like to specify a particular correlation

[R] random effects structure in lme

Dear all, In using the lme() function from the nlme package, I would like to specify a particular correlation structure for the random effects. For instance, for a 3 by 3 matrix, say, I am interested in assuming a `quasi-diagonal' matrix having zero entries in all but the first column (and the

R: [R] adding trend to an arima model

As arima.sim() simulates stationary ARMA errors if your underlying model is additive I think you can type, for instance, just: x-1:100 #time variable mu-10+.5*x #linear trend y-arima.sim(length(x), model=list(ar=.5, ma=-.3),sd=25)+mu arima(y, order=c(1,0,1),include.mean=TRUE,xreg=x) best, vito

R: [R] Change the result data

as.vector is a possible, simple solution Also use rep() on dimnames(hec.data)[[1]] to get the names vector with correct length a-matrix(1:15,ncol=5) a [,1] [,2] [,3] [,4] [,5] [1,]147 10 13 [2,]258 11 14 [3,]369 12 15 as.vector(a) [1] 1 2

R: [R] `bivariate apply' - Summary

(x,y) f( unlist(x), unlist(y) ) matrix( mapply( ff, rep(x,rep(k,k)), rep(x,k) ), k, k ) } f4-function(x,FUN,...){ #author: vito muggeo require(gregmisc) a-combinations(ncol(x),2) r-list() for(i in 1:nrow(a)){r[[length(r)+1]]- x[,a[i,]]} ris-matrix(1,ncol(x),ncol(x)) ris1

[R] `bivariate apply'

dear all, Given a matrix A, say, I would like to apply a bivariate function to each combination of its colums. That is if myfun-function(x,y)cor(x,y) #computes simple correlation of two vectors x and y then the results should be something similar to cor(A). I tried with mapply, outer,...but

R: [R] nls arguments

Dear chenu, I am not going to see your code with attention (also I do not understand the `*' symbol you used), however it looks a changepoint-type problem. The package segmented (on CRAN) uses a piecewise linear parameterization to fit regression models with breakpoints. Hope this helps, best,

R: [R] merging variable and character

See ?paste plot(x[,i], ylab=paste(series,i,sep= )) does what you want. I think it should be a good idea to read (at least) some of several official manuals, contributed doc, books on R. best, vito - Original Message - From: STOLIAROFF VINCENT [EMAIL PROTECTED] To: [EMAIL PROTECTED]

R: [R] Correction for first order autocorrelation in OLS residuals

The Cochrane Orcutt is probably an outdated approach to deal with autocorrelation and it is rather easy to write code. Why don't you use a direct likelihood-based approach? For gaussian data see the arima() function in ts package, or the Jim Lindsey's packages (for instance the gar() function in

R: [R] help with lme()

Dear Bill, I am not a lme-expert, but I believe the PinheiroBates' book is rather clear here. However you know that a lme model is, for instance fixed= y~x1+x2 and random=y~x1|group and you can fit it by ML or REML. If you are interested in testing for x2 by means the LRT (namely by comparing

[R] using variables in obj$model

dear all, for some reason I am intersted in updating a glm taking variables from its model argument, namely: dati-data.frame(y=runif(10),x=1:10) obj-glm(y~x,data=dati) obj$model[,c(A,a:b)]-cbind(rnorm(10),runif(10)) names(obj$model) [1] y x A a:b update(obj,.~.+A,data=obj$model) #it

R: [R] simulate binary data from a logistic regression model

You can use the followings: lp- -5+. #linear predictor y-rbinom(length(lp), size, plogis(lp)) Note that size means your denominator in proportions to be simulated. For instance, if you want binary data use size=1. best, vito - Original Message - From: Michele Grassi [EMAIL

[R] [R-pkgs] new package: segmented

A few days ago I uploaded to CRAN a new package called segmented. The package contains functions to fit (generalized) linear models with `segmented' (or `broken-line' or `piecewise linear') relationships between the response and one or more explanatory variables according to methodology described

[R] pmatch questions

Dear all, Below there are two, simple - I suppose, questions on using pmatch(): pmatch(xx, c(cc,xxa)) [1] 2 pmatch(a, c(cc,xxa)) [1] NA pmatch(xx, c(cc,xxa,xxb)) [1] NA I would like that the second call returns also 2, and the third call returns c(2,3) is it possible? many thanks vito

R: [R] r editors

In addition to WinEdt and (X)Emacs of course, See also, http://www.crimsoneditor.com/ (freeware) http://www.jedit.org/ (freeware) http://www.editpadpro.com/ (non-freeware) I know that there exist Syntax Highlighting files for R, but I don't know where you can find them. best, vito -

[R] models with I(1) errors

Dear all, I'm interested in fitting time-series linear models with I(1) errors. Namely given y_t=a+b*t+u_t the random term u_t are such that u_t-u_{t-1}=e_t~iid N(0,\sigma) Please, could anyone suggest me any reference (book, article, R functions) dealing with such models? Many thanks in

R: [R] extracting columns with NA's

If I have understood what you mean, you can use the (surely non-optimal) code: #build a matrix A-matrix(1:20,nrow=5) A[2,4]-NA A[,3]-rep(NA,nrow(A)) #count the missing value in each column fi-apply(A,2,function(x)sum(is.na(x))) #exclude column(s) having a number of NA equal to nrow(A). Of

[R] package documentation

Dear all, I writing my first package and everything seems to work (at least up to now) However when I try to build documentation (.dvi or .pdf), using Rcmd Rd2dvi.sh --pdf mypack.Rd I get a mypack.pdf whose title is R documentation of mypack.Rd instead of The mypack package as it should be, is it

R: [R] gam and concurvity

As someone (Simon Wood, for instance) could explain much better and as it is stressed in the help files of the mgcv pakage (the package including the gam() function) gam in R is not a clone of gam in S+. S+ uses backfitting while R uses penalized splines (see the references inside gam()

R: [R] Hosmer- Lemeshow test

As far as I know, the Hosmer-Lemeshow test is not a good choice, as it depends on the groups to be formed to compute the statistic. The design and/or hmisch packages by F.Harrell should include some alternative GoF statistics, but I didn't try. Some weeks ago there has been a similar message

R: [R] Grouping binary data

Dear Henric, The following paper deals with goodness-of-fit test for sparse (and even binary) data: Kuss O. Global goodness-of-fit tests in logistic regression with sparse data, Statist Med, 2002, 21:3789-3801. It should not too hard to write code for some non-standard and (probably under-used)

Re: [R] Combining the components of a character vector

x - c(Bob, loves, Sally) paste(x,collapse= ) [1] Bob loves Sally best, vito - Original Message - From: John Miyamoto [EMAIL PROTECTED] To: R discussion group [EMAIL PROTECTED] Sent: Thursday, April 03, 2003 1:54 AM Subject: [R] Combining the components of a character vector Dear

[R] simulating 'non-standard' survival data

Dear all, I'm looking for someone that help me to write an R function to simulate survival data under complex situations, namely time-varying hazard ratio, marginal distribution of survival times and covariates. The algorithm is described in the reference below and it should be not very difficult

[R] about the mini-distribution (R for Win )

Dear all, I've just downloaded the miniR.exe, miniR-1.bin,,miniR-8.bin in order to update R to 1.6.2 (I cannot download the file of 19Mb !). The installation seems to have worked fine, but I can't find in the /bin dir several files that were, for instance, in the /bin directory in the version

[R] error using try() and coxph()

Dear all, I'm experiencing some problems in using the following code to perform simulations with Cox models (using coxph() in survival package): for(i in 1:1000){ dd-#simulate the dataframe o0-try(coxph(Surv(start,stop,status)~x,data=dd))

[R] updating glm

Dear all, my function fn() (see code below) just takes a glm object and updates it by including a function of a specified variable in dataframe. x-1:50 y-rnorm(50) d-data.frame(yy=y,xx=x);rm(x,y) o-glm(yy~xx,data=d) fn(obj=o,x=xx) Call: glm(formula = yy ~ xx + x1, data = obj$data)

Re: [R] testing slope

Hi, in order to avoid the parameter of x:v3 to be estimated, you can re-formulate the model. You fitted y~1+x+v2+v3+x:v2+x:v3 with 6 df Now you would fit y~0+v1+v2+v3+x:v1+x:v2) with 5 df as x:v3 is not included in the model, i.e. the parameter x:v3 is constrained to be zero. best, vito A toy

Re: (v2) [R] quadratic trends and changes in slopes (R-help digest, Vol 1 #52 - 16 msgs)

It is well-known that change-point estimation is a non-trivial task. You could find interesting the followings Pastor and Guallar 1998. use of two-segmented logistic regression to estimate changepoint in epidemiological studies Am J Epid, 148, 631-642 Goetghebeur and Pocock 1995 detection and