[R] matrix^(-1/2)

2009-12-04 Thread Greg Dropkin
re [R] matrix^(-1/2) re the discussion in November on this thread. I don't know about expm but the problem must be equivalent to solve(B^(1/2)) and a solution will exist iff B is invertible and has a square root A with A%*%A = B. For 2x2 matrices necessary and sufficient conditions for B to have a

[R] scaled Schoenfeld residuals

2009-09-23 Thread Greg Dropkin
reason for the offset #or am I missing something? thanks Greg Greg Dropkin gr...@gn.apc.org __ 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.ht

Re: [R] scaled Schoenfeld residuals

2009-09-24 Thread Greg Dropkin
] [ author ] On Wed, 23 Sep 2009, Greg Dropkin wrote: > hi > > sorry if this has been discussed before, but I'm wondering why the scaled > Schoenfeld residuals do not follow the defining formula for obtaining them > from the ordinary Schoenfeld residuals, but are instead offse

[R] score test statistic in logistic regression

2009-03-03 Thread Greg Dropkin
re post from bkelcey at umich.edu bkelcey at umich.edu Wed Feb 27 15:09:48 CET 2008 If the response y is given as the proportion of successes out of n trials, and y, n, p, x, and z are vectors of length M, and the model is logit(p) = b0 + b1*x + b2*z then for the score test for the null hypothesi

[R] bsts posterior distributions

2019-01-11 Thread Greg Dropkin
hello I couldn't find previous posts to answer this, but please point me to any. I am trying to understand bsts, starting with no regressors Here is code which appears to mimic bsts, producing graphs similar to the model plot, but gives a rather different posterior distribution for the parameter

[R] gam plots and seWithMean

2010-10-21 Thread Greg Dropkin
hello I'm learning mgcv and would like to obtain numerical output corresponding to plot.gam. I can do so when seWithMean=FALSE (the default) but only approximately when seWithMean=TRUE. Can anyone show how to obtain the exact values? Alternatively, can you clarify the explanation in the manual

Re: [R] gam plots and seWithMean

2010-10-25 Thread Greg Dropkin
hadn't realised the answer would be in the source code! anyway, this appears to work. The only difference is in the last section. greg -- library(mgcv) #simulate some data x1<-runif(500) x2<-rnorm(500) x3<-rpois(500,3) d<-runif(500) t<-runif(500,20,50) linp<--6.5+x1+2*x2-x3+2*exp(-2*d)*sin(2*p

[R] Cubic splines in package "mgcv"

2011-08-16 Thread Greg Dropkin
re: Cubic splines in package "mgcv" I don't have access to Gu (2002) but clearly the function R(x,z) defined on p126 of Simon Wood's book is piecewise quartic, not piecewise cubic. Like Kunio Takezawa (below) I was puzzled by the word "cubic" on p126. As Simon Wood writes, this basis is not actu

[R] mgcv summary.gam

2013-06-21 Thread Greg Dropkin
hi I'm trying to understand (a little) the code behind summary.gam, and have the Biometrika article referred to in the Help. But am stuck early on. In the code, starting at line 167: if (est.disp) rdf <- residual.df else rdf <- -1 res <- testStat(p, Xt, V, df[i], type = p.type, res.df = rdf)

Re: [R] gam and optim

2013-09-20 Thread Greg Dropkin
please ignore this, I see the error. greg > hi > > probably a silly mistake, but I expected gam to minimise the penalised > deviance. > > thanks > > greg > > set.seed(1) > library(mgcv) > x<-runif(100) > lp<-exp(-2*x)*sin(8*x) > y<-rpois(100,exp(lp)) > plot(x,y) > m1<-gam(y~s(x),poisson) > points

[R] gam and optim

2013-09-20 Thread Greg Dropkin
hi probably a silly mistake, but I expected gam to minimise the penalised deviance. thanks greg set.seed(1) library(mgcv) x<-runif(100) lp<-exp(-2*x)*sin(8*x) y<-rpois(100,exp(lp)) plot(x,y) m1<-gam(y~s(x),poisson) points(x,exp(lp),pch=16,col="green3") points(x,fitted(m1),pch=16,cex=0.5,col="bl

Re: [R] gam and optim

2013-09-22 Thread Greg Dropkin
just to clarify how I see the error, it was the mis-definition of the penalty term in the function dv. The following code corrects this error. What is actually being minimised at this step is the penalised deviance conditional on the smoothing parameter. A second issue is that the optim default ("N

[R] mgcv: distribution of dev with Poisson data

2014-02-04 Thread Greg Dropkin
mgcv: distribution of dev hi I can't tell if this is a simple error. I'm puzzled by the distribution of dev when fitting a gam to Poisson generated data. I expected dev to be approximately chi-squared on residual d.f., i.e. about 1000 in each case below. In particular, the low values in the 3rd a

Re: [R] mgcv: distribution of dev with Poisson data

2014-02-05 Thread Greg Dropkin
roximation breaks down. So yes, the > Pearson estimator is probably a better bet in the latter case. > > best, > Simom > > ps. Note that the chi squared approximation for the *difference* in > deviance between two nested models does not suffer from this problem. > > > O

Re: [R] mgcv: distribution of dev with Poisson data

2014-02-05 Thread Greg Dropkin
h > these sort of data, and looking back at my list of open issues I see > that `someone' was you as well! It's possible that moving to a Pearson > estimator of the scale will solve that problem too. > > Thanks for this... very helpful. > > best, > Simon > >

Re: [R] mgcv: distribution of dev with Poisson data

2014-02-05 Thread Greg Dropkin
p172, another reason why I was puzzled that the default estimate used dev instead. thanks Greg > On 05/02/2014 12:56, Greg Dropkin wrote: >> thanks Simon >> also, it appears at least with ML that the default scale estimate with quasipoisson (i.e. using dev) is the scale which minim

Re: [R] REML - quasipoisson

2012-10-01 Thread Greg Dropkin
=TRUE)$val > ldsq<-sum(log(eigSq[1:16])) > F3q<-(ldhsq-ldsq)/2 > > #4th term > Mp=3 > F4q<-Mp/2*log(2*pi*phiq) > > F1q+F3q-F4q > m3$gcv > > #quite different > > #but if phiq is replaced by the Pearson est

[R] REML - quasipoisson

2012-09-25 Thread Greg Dropkin
sum(m3$edf)) F1p<-Dpq/(2*phip) F3p<-F3q #third term independent of scale F4p<-Mp/2*log(2*pi*phip) F1p+F3p-F4p m3$gcv #closer but still wrong #is there a value of phi which makes this work? f<-function(t) (Dpq/(2*t) + F3q + Mp/2*log(2*pi*t) - m3$gcv)^2 optimise(f,interval=c(0.5,1.5))$min-&

[R] mgcv: REML increase with new predictor

2012-02-24 Thread Greg Dropkin
hi Simon this follows on from the example where gcv increased unexpectedly with increasing basis dimension. I'm now looking at t2 tensor splines with REML, and find that the REML score can increase when adding a new predictor. Again, this seems odd. thanks Greg library(mgcv) set.seed(0) #simula

[R] mgcv: increasing basis dimension

2012-02-13 Thread Greg Dropkin
hi Using a ts or tprs basis, I expected gcv to decrease when increasing the basis dimension, as I thought this would minimise gcv over a larger subspace. But gcv increased. Here's an example. thanks for any comments. greg #simulate some data set.seed(0) x1<-runif(500) x2<-rnorm(500) x3<-rpois(50

Re: [R] mgcv: increasing basis dimension

2012-02-14 Thread Greg Dropkin
thanks Simon I'll upgrade R to try t2. The data I'm actually analysing requires scaled Poisson so I don't think REML is an option. thanks Greg On 14/02/12 11:22 Simon Wood wrote: That's interesting. Playing with the example, it doesn't seem to be a local minimum. I think that this happens beca