Re: [R] RODBC and sqlColumns

[Prof Brian Ripley]

AFAIK "." is not a valid part of an SQL table name. I think the help files 
are perfectly clear as to what is supported:

  sqtable: character: a database table name accessible from the
   connected dsn.

Why do you think "X.test" is a `database table name'?

On Tue, 16 Aug 2005, Ben Stabler wrote:

> I have a Postgres database that I am connecting to with the Postgres
> ODBC driver on Windows XP in R 2.1.0.  In the database is a database
> with two schemas (public and X).  With RODBC (1.1-4) , I can connect to
> the database and get the tables with sqlTables(db).  I can query tables
> in the schema with sqlQuery("SELECT * FROM X.test").  However, I can't
> get the columns in table X.test with sqlColumns(db,"X.test") //it
> returns
>
> Error in sqlColumns(db, "X.test") : 'X.test': table not found on channel
>
> If I do
>
> sqlColumns(db, "test") it returns
> [1] TABLE_QUALIFIER   TABLE_OWNER   TABLE_NAMECOLUMN_NAME
> DATA_TYPE
> [6] TYPE_NAME PRECISION LENGTHSCALE
> RADIX
> [11] NULLABLE  REMARKS   COLUMN_DEFSQL_DATA_TYPE
> SQL_DATETIME_SUB
> [16] CHAR_OCTET_LENGTH ORDINAL_POSITION  IS_NULLABLE   DISPLAY_SIZE
> FIELD_TYPE
> <0 rows> (or 0-length row.names)
>
> But there is no test table defined anywhere else but the X schema.  If I
> do sqlSave(db,aDataFrame,"X.test",T,F), it says test already defined. If
> I change the aDataFrame to be different than the fields actually in the
> data, then R starts to create a new table but returns
>
> Error in sqlColumns(db, "X.test") : 'X.test': table not found on channel
>
> It seems to be having problems with what is returned by the
> columns.since
>
> Error in sqlSave(db, aDataFrame, "X.test", T, F) :
>[RODBC] ERROR: Could not SQLExecDirectS1000 7 ERROR:  relation
> "test" already exists
>
> but if I change the input table to be differentthen R can create the
> table, but fails to populate it.  I checked the db in PgAdmin and the
> table is created by the sqlSave call.  All this stuff works if I don't
> use a schema "schema.table".  So it appears there is something wrong in
> some place dealing with understanding the columns for tables in schemas.
>
> Any ideas?  Any help would be much appreciated.  Thank you.

-- 
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] two-level poisson, again

[Henric Nilsson]

On On, 2005-08-17, 08:52, Shige Song skrev:
> Hi,
>
> I compare results of a simple two-level poisson estimated using lmer
> and those estimated using MLwiN and Stata (v.9).
>
> In R, I trype:
> ---
> m2 <- lmer(.D ~ offset(log(.Y)) + (1|pcid2) + educy + agri, male, poisson)
> ---
>
> In Stata, I type:
> ---
> xtpois _D educy agri, e(_Y) i(pcid2)

You're not fitting the same model! `lmer' uses Gaussian random effects and
the default for `xtpois' is gamma random effects.

Also, note that even if you'd specified a Gaussian random effect (through
a `normal' to the right of the `,' in your `xtpois' call) the same fitting
criterion is not used since `xtpois' uses adaptive Gauss-Hermite
quadrature and `lmer' defaults to PQL.

For comparable results, try the following:

m2 <- lmer(.D ~ offset(log(.Y)) + (1|pcid2) + educy + agri, male, poisson,
method = "AGQ")

xtpois _D educy agri, e(_Y) i(pcid2) re normal

You may also want to try Göran Broström's `glmmML' package.


HTH,
Henric

>
> Results using R look like:
> ---
> ..
> Random effects:
>  GroupsNameVarianceStd.Dev.
>   pcid2 (Intercept)   5e-10  2.2361e-05
> # of obs: 25360, groups: pcid2, 174
>
> Estimated scale (compare to 1)  7.190793
>
> Fixed effects:
>Estimate  Std. Error z value  Pr(>|z|)
> (Intercept) -3.28548086  0.00408771 -803.75 < 2.2e-16 ***
> educy0.00507475  0.00039616   12.81 < 2.2e-16 ***
> agri 0.01346887  0.003348374.02 5.758e-05 ***
> ..
> --
>
> Results using Stata look like:
>
> --
>   _D |  Coef.   Std. Err.  zP>|z| [95% Conf.
> Interval]
> -+
>educy |   .0120431   .000444127.12   0.000 .0111725
> .0129136
> agri |   .0293177   .0035586 8.24   0.000  .022343
> .0362924
>_cons |  -3.325073   .0076275  -435.93   0.000-3.340023
> -3.310123
>   _Y | (exposure)
> -+
> /lnalpha |  -4.982977   .1156474 -5.209641
> -4.756312
> -+
>alpha |   .0068536   .0007926  .0054636
> .0085973
> --
>
>
> As you can see, the discrepency is significant! And results using
> MLwiN agree with Stata. Any help will be greately appreciated!
>
> Shige
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] Mixed Effects Model Power Calculations

[Henric Nilsson]

On Ti, 2005-08-16, 21:17, [EMAIL PROTECTED] skrev:

> Is there an R package available that would facilitate doing a power/sample
> size analysis for linear mixed effects models?

I'm not aware of such a package (others might be...).

When it comes to sample size calculations, especially for tricky designs
and/or advanced methodology, simulation is usually the best approach. An
example using `lme' can be found at

http://maven.smith.edu/~nhorton/R/


HTH,
Henric

> I have seen the Java applets made available by Russell Length which would
> seem to be able to handle most any lme, but there is little documentation
> and it's not clear how the models need to be formulated.
>
> Rick B.
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] Mixed Effects Model Power Calculations

[Dimitris Rizopoulos]

I don't know what specific application Rick has in mind, but if there 
is possibility of missing values (which is common, e.g., in 
longitudinal studies) then this should also be taken into account in 
the power calculations.

Best,
Dimitris


Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/16/336899
Fax: +32/16/337015
Web: http://www.med.kuleuven.be/biostat/
 http://www.student.kuleuven.be/~m0390867/dimitris.htm


- Original Message - 
From: "Henric Nilsson" <[EMAIL PROTECTED]>
To: <[EMAIL PROTECTED]>
Cc: 
Sent: Wednesday, August 17, 2005 10:42 AM
Subject: Re: [R] Mixed Effects Model Power Calculations


>
> On Ti, 2005-08-16, 21:17, [EMAIL PROTECTED] skrev:
>
>> Is there an R package available that would facilitate doing a 
>> power/sample
>> size analysis for linear mixed effects models?
>
> I'm not aware of such a package (others might be...).
>
> When it comes to sample size calculations, especially for tricky 
> designs
> and/or advanced methodology, simulation is usually the best 
> approach. An
> example using `lme' can be found at
>
> http://maven.smith.edu/~nhorton/R/
>
>
> HTH,
> Henric
>
>> I have seen the Java applets made available by Russell Length which 
>> would
>> seem to be able to handle most any lme, but there is little 
>> documentation
>> and it's not clear how the models need to be formulated.
>>
>> Rick B.
>>
>> __
>> R-help@stat.math.ethz.ch mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide!
>> http://www.R-project.org/posting-guide.html
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] Two-level Poisson model with cross classified random effects

[Renaud Lancelot]

Shige Song a écrit :
> Dear All,
> 
> I have two-level data with individual as level-1, birth cohort and
> community as level-2. All the level-2 covariates are generated from
> the level-1 covariates by cross-classifying by cohort and community.
> 
>>From what I read, an ordinary three-level model with individual nesed
> within birth cohort nested within community, or individual nested
> within community nested within birth cohort do not work well, neither
> do model with individual nested within community by cohort. The right
> way to go is to estimate a two-level model with two separate random
> effects: within cohort and within community. The question I want to
> ask is: how to do this using lmer?
> 
> I tried the following for a simple unconditioal model:
> 
> m1 <- lmer(count ~ offset(log(total)) + (1|comm) + (1|cohort), data, poisson)
> 
> where "count" is the dependent variable, "total" is the exposure
> variable, "comm" is the community ID, and "cohort" is the birth cohort
> ID. Will this be suffice? I got really smalle randome intercept
> (5.e-10 for community and 4.4226e-05 for cohort), which got me a
> bit nervous.
> 
> Thanks!
> 
> Best,
> Shige
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
> 

The default method to fit generalized mixed-effect models with lme4 is 
"PQL" which is fast but not very accurate for the random effects (they 
might be underestimated). Try other methods ("Laplace" or "AGQ") to see 
if you get different results.

Best,

Renaud

-- 
Dr Renaud Lancelot, vétérinaire
Projet FSP régional épidémiologie vétérinaire
C/0 Ambassade de France - SCAC
BP 834 Antananarivo 101 - Madagascar

e-mail: [EMAIL PROTECTED]
tel.:   +261 32 40 165 53 (cell)
 +261 20 22 665 36 ext. 225 (work)

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] Mixed Effects Model Power Calculations

[Shige Song]

Hi Dimitris,

Thank you so much! That really helps!

Shige

On 8/17/05, Dimitris Rizopoulos <[EMAIL PROTECTED]> wrote:
> I don't know what specific application Rick has in mind, but if there
> is possibility of missing values (which is common, e.g., in
> longitudinal studies) then this should also be taken into account in
> the power calculations.
> 
> Best,
> Dimitris
> 
> 
> Dimitris Rizopoulos
> Ph.D. Student
> Biostatistical Centre
> School of Public Health
> Catholic University of Leuven
> 
> Address: Kapucijnenvoer 35, Leuven, Belgium
> Tel: +32/16/336899
> Fax: +32/16/337015
> Web: http://www.med.kuleuven.be/biostat/
>  http://www.student.kuleuven.be/~m0390867/dimitris.htm
> 
> 
> - Original Message -
> From: "Henric Nilsson" <[EMAIL PROTECTED]>
> To: <[EMAIL PROTECTED]>
> Cc: 
> Sent: Wednesday, August 17, 2005 10:42 AM
> Subject: Re: [R] Mixed Effects Model Power Calculations
> 
> 
> >
> > On Ti, 2005-08-16, 21:17, [EMAIL PROTECTED] skrev:
> >
> >> Is there an R package available that would facilitate doing a
> >> power/sample
> >> size analysis for linear mixed effects models?
> >
> > I'm not aware of such a package (others might be...).
> >
> > When it comes to sample size calculations, especially for tricky
> > designs
> > and/or advanced methodology, simulation is usually the best
> > approach. An
> > example using `lme' can be found at
> >
> > http://maven.smith.edu/~nhorton/R/
> >
> >
> > HTH,
> > Henric
> >
> >> I have seen the Java applets made available by Russell Length which
> >> would
> >> seem to be able to handle most any lme, but there is little
> >> documentation
> >> and it's not clear how the models need to be formulated.
> >>
> >> Rick B.
> >>
> >> __
> >> R-help@stat.math.ethz.ch mailing list
> >> https://stat.ethz.ch/mailman/listinfo/r-help
> >> PLEASE do read the posting guide!
> >> http://www.R-project.org/posting-guide.html
> >
> > __
> > R-help@stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide!
> > http://www.R-project.org/posting-guide.html
> >
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] Two-level Poisson model with cross classified random effects

[Shige Song]

Dear Renaud,

Thank you very much! So now lme4 has both laplace and adaptive
quadrature? Wow, that's really impressive!

Shige

On 8/17/05, Renaud Lancelot <[EMAIL PROTECTED]> wrote:
> Shige Song a écrit :
> > Dear All,
> >
> > I have two-level data with individual as level-1, birth cohort and
> > community as level-2. All the level-2 covariates are generated from
> > the level-1 covariates by cross-classifying by cohort and community.
> >
> >>From what I read, an ordinary three-level model with individual nesed
> > within birth cohort nested within community, or individual nested
> > within community nested within birth cohort do not work well, neither
> > do model with individual nested within community by cohort. The right
> > way to go is to estimate a two-level model with two separate random
> > effects: within cohort and within community. The question I want to
> > ask is: how to do this using lmer?
> >
> > I tried the following for a simple unconditioal model:
> >
> > m1 <- lmer(count ~ offset(log(total)) + (1|comm) + (1|cohort), data, 
> > poisson)
> >
> > where "count" is the dependent variable, "total" is the exposure
> > variable, "comm" is the community ID, and "cohort" is the birth cohort
> > ID. Will this be suffice? I got really smalle randome intercept
> > (5.e-10 for community and 4.4226e-05 for cohort), which got me a
> > bit nervous.
> >
> > Thanks!
> >
> > Best,
> > Shige
> >
> > __
> > R-help@stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide! 
> > http://www.R-project.org/posting-guide.html
> >
> 
> The default method to fit generalized mixed-effect models with lme4 is
> "PQL" which is fast but not very accurate for the random effects (they
> might be underestimated). Try other methods ("Laplace" or "AGQ") to see
> if you get different results.
> 
> Best,
> 
> Renaud
> 
> --
> Dr Renaud Lancelot, vétérinaire
> Projet FSP régional épidémiologie vétérinaire
> C/0 Ambassade de France - SCAC
> BP 834 Antananarivo 101 - Madagascar
> 
> e-mail: [EMAIL PROTECTED]
> tel.:   +261 32 40 165 53 (cell)
>  +261 20 22 665 36 ext. 225 (work)
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] power of a matrix

[Rau, Roland]

Dear all,

I have a population with three age-classes, at time t=0 the population
is:
n.zero <- c(1,0,0)
I have a transition matrix A which denotes "fertility" and "survival":

A <- matrix(c(0,1,5, 0.3,0,0, 0,0.5,0), ncol=3, byrow=TRUE)

To obtain the population at t=1, I calculate:
A %*% n.zero

To obtain the population t=2, I calculate:
A %*% (A %*% n.zero)

... and so on ...

I thought now to obtain the population at time x, I should simply do:
A^x %*% n.zero

But this, of course, does not work since the following two statements
are not equivalent for matrices:
A %*% A
A * A

Is there something like a "powermatrix"-function?

Thanks,
Roland

P.S. The example is taken from:
Example 3.1 ("A linear, time-invariant model") from Keyfitz/Caswell:
"Applied Mathematical Demography", 3rd Edition, 2005, p. 52f


+
This mail has been sent through the MPI for Demographic Rese...{{dropped}}

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] How to assess significance of random effect in lme4

[Shige Song]

Dear All,

With kind help from several friends on the list, I am getting close.
Now here are something interesting I just realized: for random
effects, lmer reports standard deviation instead of standard error! Is
there a hidden option that tells lmer to report standard error of
random effects, like most other multilevel or mixed modeling software,
so that we can say something like "randome effect for xxx is
significant, while randome effect for xxx is not significant"? Thanks!

Best,
Shige

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] x-axis binning

[Samuel Kemp]

Hi,

I am wanting to plot an x,y data.frame. But, I would like to have bins 
on the x-axis. Is there an easy way of doing this?

Thanks in advance,

Sam.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] R: characters

[Clark Allan]

hi all 

assume that i have the following table

a=rbind(c(T,T,F),c(F,F,T))
> a
  [,1]  [,2]  [,3]
[1,]  TRUE  TRUE FALSE
[2,] FALSE FALSE  TRUE

I would like to change all the FALSE entries to a blank. how can i do
this?

i could simply use

a[a==F]=""
a

but then how would i remove the " " from the entries.

i know that this should be very easy!!!

/
allan__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] R: characters

[Dimitris Rizopoulos]

you could convert it to a data.frame, i.e.,

a <- rbind(c(T, T, F), c(F, F, T))
a[!a] <- ""
as.data.frame(a)


I hope it helps.

Best,
Dimitris


Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/16/336899
Fax: +32/16/337015
Web: http://www.med.kuleuven.be/biostat/
 http://www.student.kuleuven.be/~m0390867/dimitris.htm


- Original Message - 
From: "Clark Allan" <[EMAIL PROTECTED]>
To: 
Sent: Wednesday, August 17, 2005 12:35 PM
Subject: [R] R: characters


> hi all
>
> assume that i have the following table
>
> a=rbind(c(T,T,F),c(F,F,T))
>> a
>  [,1]  [,2]  [,3]
> [1,]  TRUE  TRUE FALSE
> [2,] FALSE FALSE  TRUE
>
> I would like to change all the FALSE entries to a blank. how can i 
> do
> this?
>
> i could simply use
>
> a[a==F]=""
> a
>
> but then how would i remove the " " from the entries.
>
> i know that this should be very easy!!!
>
> /
> allan





> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] x-axis binning

[Petr Pikal]

Hi

On 17 Aug 2005 at 11:34, Samuel Kemp wrote:

> Hi,
> 
> I am wanting to plot an x,y data.frame. But, I would like to have bins
> on the x-axis. Is there an easy way of doing this?

Yes.

Do you want to know how?

?cut

is probably what you want, but as you did not tell us any further 
details it is hard to decipher what you exactly want.

HTH
Petr



> 
> Thanks in advance,
> 
> Sam.
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html

Petr Pikal
[EMAIL PROTECTED]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] How to assess significance of random effect in lme4

[Doran, Harold]

You can extract the posterior variance of the random effect from the bVar slot 
of the fitted lmer model. It is not a hidden option, but a part of the fitted 
model. It just doesn't show up when you use summary().

Look at the structure of your object to see what is available using str().

However, your comment below seems to imply that it is incorrect for lmer to 
report SDs instead of the standard error, which is not true. That is a quantity 
of direct interest.

Other multilevel programs report the same exact statistics (for the most part). 
For instance, HLM reports the variances as well. If you want the posterior 
variance of an HLM model you need to extract it.


-Original Message-
From:   [EMAIL PROTECTED] on behalf of Shige Song
Sent:   Wed 8/17/2005 6:30 AM
To: r-help@stat.math.ethz.ch
Cc: 
Subject:[R] How to assess significance of random effect in lme4

Dear All,

With kind help from several friends on the list, I am getting close.
Now here are something interesting I just realized: for random
effects, lmer reports standard deviation instead of standard error! Is
there a hidden option that tells lmer to report standard error of
random effects, like most other multilevel or mixed modeling software,
so that we can say something like "randome effect for xxx is
significant, while randome effect for xxx is not significant"? Thanks!

Best,
Shige

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html




[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] power of a matrix

[Dimitris Rizopoulos]

look at function ?mtx.exp() in the Malmig package, e.g.,

A <- matrix(c(0,1,5, 0.3,0,0, 0,0.5,0), ncol=3, byrow=TRUE)
n.zero <- c(1,0,0)
##
library(Malmig)

all.equal(A %*% (A %*% n.zero), mtx.exp(A, 2) %*% n.zero)
all.equal(A %*% (A %*% (A %*% n.zero)), mtx.exp(A, 3) %*% n.zero)

mtx.exp(A, 15) %*% n.zero


I hope it helps.

Best,
Dimitris


Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/16/336899
Fax: +32/16/337015
Web: http://www.med.kuleuven.be/biostat/
 http://www.student.kuleuven.be/~m0390867/dimitris.htm


- Original Message - 
From: "Rau, Roland" <[EMAIL PROTECTED]>
To: "R-Help" 
Sent: Wednesday, August 17, 2005 12:15 PM
Subject: [R] power of a matrix


> Dear all,
>
> I have a population with three age-classes, at time t=0 the 
> population
> is:
> n.zero <- c(1,0,0)
> I have a transition matrix A which denotes "fertility" and 
> "survival":
>
> A <- matrix(c(0,1,5, 0.3,0,0, 0,0.5,0), ncol=3, byrow=TRUE)
>
> To obtain the population at t=1, I calculate:
> A %*% n.zero
>
> To obtain the population t=2, I calculate:
> A %*% (A %*% n.zero)
>
> ... and so on ...
>
> I thought now to obtain the population at time x, I should simply 
> do:
> A^x %*% n.zero
>
> But this, of course, does not work since the following two 
> statements
> are not equivalent for matrices:
> A %*% A
> A * A
>
> Is there something like a "powermatrix"-function?
>
> Thanks,
> Roland
>
> P.S. The example is taken from:
> Example 3.1 ("A linear, time-invariant model") from Keyfitz/Caswell:
> "Applied Mathematical Demography", 3rd Edition, 2005, p. 52f
>
>
> +
> This mail has been sent through the MPI for Demographic 
> Rese...{{dropped}}
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] How to assess significance of random effect in lme4

[Shige Song]

Hi Harold,

Thanks for the reply. I looked at my outputs using str() as you
suggested, here is the part you mentioned:

  ..@ bVar :List of 2
  .. ..$ commu  : num [1, 1, 1:29] 5e-10 5e-10 5e-10 5e-10 5e-10 ...
  .. ..$ bcohort: num [1, 1, 1:6] 1.05e-05 7.45e-06 6.53e-06 8.25e-06
7.11e-06 ...

where commu and bcohort are the two second-level units. Are these
standard errors? Why the second vector contains a series of different
numbers?

Thanks!

Shige

On 8/17/05, Doran, Harold <[EMAIL PROTECTED]> wrote:
>  
> 
> You can extract the posterior variance of the random effect from the bVar
> slot of the fitted lmer model. It is not a hidden option, but a part of the
> fitted model. It just doesn't show up when you use summary().
>  
>  Look at the structure of your object to see what is available using str().
>  
>  However, your comment below seems to imply that it is incorrect for lmer to
> report SDs instead of the standard error, which is not true. That is a
> quantity of direct interest.
>  
>  Other multilevel programs report the same exact statistics (for the most
> part). For instance, HLM reports the variances as well. If you want the
> posterior variance of an HLM model you need to extract it.
> 
>  
>  
>  -Original Message-
>  From:   [EMAIL PROTECTED] on behalf of
> Shige Song
>  Sent:   Wed 8/17/2005 6:30 AM
>  To: r-help@stat.math.ethz.ch
>  Cc:
>  Subject:[R] How to assess significance of random effect in lme4
>  
>  Dear All,
>  
>  With kind help from several friends on the list, I am getting close.
>  Now here are something interesting I just realized: for random
>  effects, lmer reports standard deviation instead of standard error! Is
>  there a hidden option that tells lmer to report standard error of
>  random effects, like most other multilevel or mixed modeling software,
>  so that we can say something like "randome effect for xxx is
>  significant, while randome effect for xxx is not significant"? Thanks!
>  
>  Best,
>  Shige
>  
>  __
>  R-help@stat.math.ethz.ch mailing list
>  https://stat.ethz.ch/mailman/listinfo/r-help
>  PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>  
>  
>  
> 
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] R: characters

[Martin Lam]

Hmm, it seems that the only difference is the use of
"as.data.frame(a)" instead of "a". Hence, the same
result can be done with:

a <- rbind(c(T, T, F), c(F, F, T))
a[a==F]=""
as.data.frame(a)

A possible drawback is that the mode changes from
logical to character when using a[a==F]=""

Instead you could use a[a==F]=NA, but that won't get
you blank lines.

Best,

Martin


--- Dimitris Rizopoulos
<[EMAIL PROTECTED]> wrote:

> you could convert it to a data.frame, i.e.,
> 
> a <- rbind(c(T, T, F), c(F, F, T))
> a[!a] <- ""
> as.data.frame(a)
> 
> 
> I hope it helps.
> 
> Best,
> Dimitris
> 
> 
> Dimitris Rizopoulos
> Ph.D. Student
> Biostatistical Centre
> School of Public Health
> Catholic University of Leuven
> 
> Address: Kapucijnenvoer 35, Leuven, Belgium
> Tel: +32/16/336899
> Fax: +32/16/337015
> Web: http://www.med.kuleuven.be/biostat/
> 
>
http://www.student.kuleuven.be/~m0390867/dimitris.htm
> 
> 
> - Original Message - 
> From: "Clark Allan" <[EMAIL PROTECTED]>
> To: 
> Sent: Wednesday, August 17, 2005 12:35 PM
> Subject: [R] R: characters
> 
> 
> > hi all
> >
> > assume that i have the following table
> >
> > a=rbind(c(T,T,F),c(F,F,T))
> >> a
> >  [,1]  [,2]  [,3]
> > [1,]  TRUE  TRUE FALSE
> > [2,] FALSE FALSE  TRUE
> >
> > I would like to change all the FALSE entries to a
> blank. how can i 
> > do
> > this?
> >
> > i could simply use
> >
> > a[a==F]=""
> > a
> >
> > but then how would i remove the " " from the
> entries.
> >
> > i know that this should be very easy!!!
> >
> > /
> > allan
> 
> 
>

> 
> 
> > __
> > R-help@stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide! 
> > http://www.R-project.org/posting-guide.html
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] R: characters

[Peter Dalgaard]

Martin Lam <[EMAIL PROTECTED]> writes:

> Hmm, it seems that the only difference is the use of
> "as.data.frame(a)" instead of "a". Hence, the same
> result can be done with:
> 
> a <- rbind(c(T, T, F), c(F, F, T))
> a[a==F]=""
> as.data.frame(a)

That's silly:

> a[!a]<-""
> noquote(a)
 [,1] [,2] [,3]
[1,] TRUE TRUE
[2,]   TRUE


-- 
   O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark  Ph:  (+45) 35327918
~~ - ([EMAIL PROTECTED])  FAX: (+45) 35327907

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] power of a matrix

[Rau, Roland]

Thank you very much! Thanks also to the authors of this function,
Vincente Canto Cassola and Martin Maechler!

This is exactly what I hoped for.

Best,
Roland



> -Original Message-
> From: Dimitris Rizopoulos 
> [mailto:[EMAIL PROTECTED] 
> Sent: Wednesday, August 17, 2005 1:38 PM
> To: Rau, Roland
> Cc: r-help@stat.math.ethz.ch
> Subject: Re: [R] power of a matrix
> 
> look at function ?mtx.exp() in the Malmig package, e.g.,
> 
> A <- matrix(c(0,1,5, 0.3,0,0, 0,0.5,0), ncol=3, byrow=TRUE)
> n.zero <- c(1,0,0)
> ##
> library(Malmig)
> 
> all.equal(A %*% (A %*% n.zero), mtx.exp(A, 2) %*% n.zero)
> all.equal(A %*% (A %*% (A %*% n.zero)), mtx.exp(A, 3) %*% n.zero)
> 
> mtx.exp(A, 15) %*% n.zero
> 
> 
> I hope it helps.
> 
> Best,
> Dimitris
> 
> 
> Dimitris Rizopoulos
> Ph.D. Student
> Biostatistical Centre
> School of Public Health
> Catholic University of Leuven
> 
> Address: Kapucijnenvoer 35, Leuven, Belgium
> Tel: +32/16/336899
> Fax: +32/16/337015
> Web: http://www.med.kuleuven.be/biostat/
>  http://www.student.kuleuven.be/~m0390867/dimitris.htm
> 
> 
> - Original Message - 
> From: "Rau, Roland" <[EMAIL PROTECTED]>
> To: "R-Help" 
> Sent: Wednesday, August 17, 2005 12:15 PM
> Subject: [R] power of a matrix
> 
> 
> > Dear all,
> >
> > I have a population with three age-classes, at time t=0 the 
> > population
> > is:
> > n.zero <- c(1,0,0)
> > I have a transition matrix A which denotes "fertility" and 
> > "survival":
> >
> > A <- matrix(c(0,1,5, 0.3,0,0, 0,0.5,0), ncol=3, byrow=TRUE)
> >
> > To obtain the population at t=1, I calculate:
> > A %*% n.zero
> >
> > To obtain the population t=2, I calculate:
> > A %*% (A %*% n.zero)
> >
> > ... and so on ...
> >
> > I thought now to obtain the population at time x, I should simply 
> > do:
> > A^x %*% n.zero
> >
> > But this, of course, does not work since the following two 
> > statements
> > are not equivalent for matrices:
> > A %*% A
> > A * A
> >
> > Is there something like a "powermatrix"-function?
> >
> > Thanks,
> > Roland
> >
> > P.S. The example is taken from:
> > Example 3.1 ("A linear, time-invariant model") from Keyfitz/Caswell:
> > "Applied Mathematical Demography", 3rd Edition, 2005, p. 52f
> >
> >
> > +
> > This mail has been sent through the MPI for Demographic 
> > Rese...{{dropped}}
> >
> > __
> > R-help@stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide! 
> > http://www.R-project.org/posting-guide.html
> > 
> 
> 


+
This mail has been sent through the MPI for Demographic Rese...{{dropped}}

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] power of a matrix

[Peter Dalgaard]

"Rau, Roland" <[EMAIL PROTECTED]> writes:

> Thank you very much! Thanks also to the authors of this function,
> Vincente Canto Cassola and Martin Maechler!
> 
> This is exactly what I hoped for.

> > look at function ?mtx.exp() in the Malmig package, e.g.,

Also, there is mexp() in the Matrix package. I'm not sure about the
relative merits. mexp() is one of the less dubious implementations of
matrix exponentials, but it does require to and from class "Matrix".
mtx.exp is a bit unfortunately named as it appears to calculate matrix
*powers* (which in this case is what you need).

-- 
   O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark  Ph:  (+45) 35327918
~~ - ([EMAIL PROTECTED])  FAX: (+45) 35327907

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] error in predict glm (new levels cause problems)

[Spencer Graves]

  Have you considered replacing the offending factor with explicit 
coding of your own choosing?  See "?contr.helmert" and in library(MASS) 
"contr.sdif" plus the "contrasts" attribute in "options", obtained, 
e.g., via 'options("contrasts")'.  Every k-level factor is by default 
converted into a set of (k-1) columns of numeric codes.  The exact 
numbers do not matter to p values obtained from anova, though they will 
matter to the coefficients estimated.

  Have you tried something like the following with "glm.nb":

 > set.seed(1)
 > DF <- data.frame(a=sample(letters[1:3], 30, replace=TRUE),
+  b=sample(LETTERS[1:3], 30, replace=TRUE),
+  y=rnorm(30))
 > with(DF, table(a, b))
b
a   A B C
   a 3 3 3
   b 2 6 3
   c 2 4 4
 > fit1 <- lm(y~a+b, DF[DF$a!="a",])
 > fit1$contrasts
$a
[1] "contr.treatment"

$b
[1] "contr.treatment"

 > options(contrasts=c(unordered="contr.helmert", ordered="contr.poly"))
 > fit2 <- lm(y~a+b, DF[DF$a!="a",])
 > fit2$contrasts
$a
[1] "contr.helmert"

$b
[1] "contr.helmert"

 > anova(fit1)
Analysis of Variance Table

Response: y
   Df Sum Sq Mean Sq F value Pr(>F)
a  1 0.0008  0.0008  0.0015 0.9698
b  2 0.6186  0.3093  0.5775 0.5719
Residuals 17 9.1050  0.5356
 > anova(fit2)
Analysis of Variance Table

Response: y
   Df Sum Sq Mean Sq F value Pr(>F)
a  1 0.0008  0.0008  0.0015 0.9698
b  2 0.6186  0.3093  0.5775 0.5719
Residuals 17 9.1050  0.5356
 >
  spencer graves

K. Steinmann wrote:

> Dear R-helpers,
> 
> I try to perform glm's with negative binomial distributed data.
> So I use the MASS library and the commands:
> model_1 = glm.nb(response ~ y1 + y2 + ...+ yi, data = data.frame)
> and
> predict(model_1, newdata = data.frame)
> 
> 
> So far, I think everything should be ok.
> 
> But when I want to perform a glm with a subset of the data,
> I run into an error message as soon as I want to predict values, based on the
> new model. The problem seems to be the reduced number of levels of one of the
> factors yi ( a categorical factor) in the subset of the original data set.
> 
> On cran search I found some related hint, that the line "mf$drop.unused.levels
> <- TRUE " in the glm (or glm.nb) function could cause the problem.
> 
> Therefore I changed the line to "mf$drop.unused.levels <- FALSE ".
> Indeed the error message disappears and when I compare the prediction of 
> model_1
> with the prediction of the model, carried out with the full data set but with
> the changed glm.nb function, I get the same predicted numbers.
> 
> However, the change of glm.nb function was more of an intuitive action, and
> since I still consider myself as a beginner of R, I don't feel comfortable.
> 
> So my questions:
> 1. Is there an easier way to solve my problem?
> 2. Do I affect the glm.nb function seriously, by changing the line mentioned
> above?
> 
> 
> Thank you for your help,
> Katharina
> 
> PS: I am working with R 2.0.0
> PPS: Concrete error message:
> "Error in model.frame.default(Terms, newdata, na.action = na.action, xlev =
> object$xlevels) :
> factor I(kanton) has new level(s) GE"
> 
> 
> 
> 
> --
> K. Steinmann
> Botanisches Institut
> Universität Basel
> CH-4056 Basel
> Switzerland
> Tel  0041 61 267 35 02
> E-mail: [EMAIL PROTECTED]
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

-- 
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

[EMAIL PROTECTED]
www.pdf.com 
Tel:  408-938-4420
Fax: 408-280-7915

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] How to assess significance of random effect in lme4

[Doran, Harold]

These are the posterior variances of the random effects (I think more
properly termed "empirical" posteriors).  Your model apparently includes
three levels of random variation (commu, bcohort, residual). The first
are the variances associated with your commu random effect and the
second are the variances associated with the bcohort random effect.

Accessing either one would require

[EMAIL PROTECTED] or [EMAIL PROTECTED]

Obviously, replace "fm" with the name of your fitted model.

-Original Message-
From: [EMAIL PROTECTED]
[mailto:[EMAIL PROTECTED] On Behalf Of Shige Song
Sent: Wednesday, August 17, 2005 7:50 AM
To: r-help@stat.math.ethz.ch
Subject: Re: [R] How to assess significance of random effect in lme4

Hi Harold,

Thanks for the reply. I looked at my outputs using str() as you
suggested, here is the part you mentioned:

  ..@ bVar :List of 2
  .. ..$ commu  : num [1, 1, 1:29] 5e-10 5e-10 5e-10 5e-10 5e-10 ...
  .. ..$ bcohort: num [1, 1, 1:6] 1.05e-05 7.45e-06 6.53e-06 8.25e-06
7.11e-06 ...

where commu and bcohort are the two second-level units. Are these
standard errors? Why the second vector contains a series of different
numbers?

Thanks!

Shige

On 8/17/05, Doran, Harold <[EMAIL PROTECTED]> wrote:
>  
> 
> You can extract the posterior variance of the random effect from the 
> bVar slot of the fitted lmer model. It is not a hidden option, but a 
> part of the fitted model. It just doesn't show up when you use
summary().
>  
>  Look at the structure of your object to see what is available using
str().
>  
>  However, your comment below seems to imply that it is incorrect for 
> lmer to report SDs instead of the standard error, which is not true. 
> That is a quantity of direct interest.
>  
>  Other multilevel programs report the same exact statistics (for the 
> most part). For instance, HLM reports the variances as well. If you 
> want the posterior variance of an HLM model you need to extract it.
> 
>  
>  
>  -Original Message-
>  From:   [EMAIL PROTECTED] on behalf of
> Shige Song
>  Sent:   Wed 8/17/2005 6:30 AM
>  To: r-help@stat.math.ethz.ch
>  Cc:
>  Subject:[R] How to assess significance of random effect in
lme4
>  
>  Dear All,
>  
>  With kind help from several friends on the list, I am getting close.
>  Now here are something interesting I just realized: for random  
> effects, lmer reports standard deviation instead of standard error! Is

> there a hidden option that tells lmer to report standard error of  
> random effects, like most other multilevel or mixed modeling software,

> so that we can say something like "randome effect for xxx is  
> significant, while randome effect for xxx is not significant"? Thanks!
>  
>  Best,
>  Shige
>  
>  __
>  R-help@stat.math.ethz.ch mailing list  
> https://stat.ethz.ch/mailman/listinfo/r-help
>  PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>  
>  
>  
> 
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] resampling Question

[Achaz von Hardenberg]

hi,
sorry for a possibly naive question but I am a bit of a beginner in R 
programming...

I wrote a function  which simulates some data and performs two different kinds 
of analyses on it. as an output I get the statistics for the two analyses 
(t-values). Now I would like to have an other function which reruns my first 
function say a 1000 times and attaches the resulting statistics in a data.frame 
so that at the end it contains 1000 rows with two columns for the two 
statistics I calculated with my first function.

I hope I was clear enugh and would be glad to anybody who can help me out!

ciao,
 
Achaz von Hardenberg
--
Centro Studi Fauna Alpina - Alpine Wildlife Research Centre
Parco Nazionale Gran Paradiso, via della Rocca 47, 10123 Torino, Italy

e-mail: [EMAIL PROTECTED]
Tel. (office): +39.011.8606212
Tel. (mobile): +39.328.8736291
Fax: +39.011.8121305
--
Open access to all papers published in the Journal of Mountain Ecology:
http://www.mountainecology.org

GSE-AIESG (Gruppo Stambecco Europa - Alpine Ibex European Specialist Group):
http://www.gseonline.org

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] resampling Question

[Peter Dalgaard]

"Achaz von Hardenberg" <[EMAIL PROTECTED]> writes:

> hi,
> sorry for a possibly naive question but I am a bit of a beginner in R 
> programming...
> 
> I wrote a function which simulates some data and performs two
> different kinds of analyses on it. as an output I get the statistics
> for the two analyses (t-values). Now I would like to have an other
> function which reruns my first function say a 1000 times and
> attaches the resulting statistics in a data.frame so that at the end
> it contains 1000 rows with two columns for the two statistics I
> calculated with my first function.
> 
> I hope I was clear enugh and would be glad to anybody who can help me out!

replicate is a good start:

> replicate(10,t.test(rexp(10),mu=1)[c("statistic","p.value")])
  [,1]  [,2]  [,3]   [,4]  [,5]   [,6]
statistic -1.552478 -1.09727  -2.053807  -1.671855 -0.1169578 0.1005037
p.value   0.1549645 0.3010120 0.07018008 0.1288855 0.9094619  0.9221477
  [,7]   [,8]  [,9]  [,10]
statistic -0.1711356 -0.933484 0.3169710 -1.136498
p.value   0.867903   0.3749356 0.758495  0.2851029

As you see, the result is a 2xn matrix. If you really need a data
frame, just use as.data.frame(t()),  or (bypassing the matrix
entirely):

> do.call("rbind",replicate(10,as.data.frame(
   t.test(rexp(10),mu=1)[c("statistic","p.value")]), simplify=FALSE))
 statistic   p.value
t   0.32566430 0.7521223
t1 -1.22741479 0.2508023
t2 -1.66792987 0.1296757
t3  1.56440274 0.1521619
t4  0.63778111 0.5395015
t5 -1.03826715 0.3262346
t6  0.09337127 0.9276542
t7  0.90166085 0.3907282
t8 -0.78164107 0.4544958
t9 -0.39766367 0.7001452


-- 
   O__   Peter Dalgaard Øster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark  Ph:  (+45) 35327918
~~ - ([EMAIL PROTECTED])  FAX: (+45) 35327907

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] Fwd: Documenting data sets with many variables

[Arne Henningsen]

On Tuesday 16 August 2005 17:26, Gavin Simpson wrote:
> On Tue, 2005-08-16 at 17:11 +0200, Arne Henningsen wrote:
> > On Tuesday 16 August 2005 14:49, Roger D. Peng wrote:
> > > Have you tried using 'promptData()' on the data frame and then
> > > just using the resulting documentation file?
> >
> > Thank you, Roger, for bringing 'promptData()' to my mind. This is really
> > a useful tool. However, in my special case my aim is to reduce the extent
> > and increase the comprehensibility of the documentation rather than to
> > reduce my effort to write the documentation.
> >
> > Any further hints are welcome!
> >
> > Thanks,
> > Arne
>
> Would it not be expedient then to ignore the \format{} section and just
> provide the information on the variables say in the \description{},
> e.g.:

That's a great idea - and so simple!
This perfectly solves my problem.
Thanks,
Arne

> This example taken from package vegan describing 2 data.frames with 44
> and 14 columns. Admittedly, none of the variables in the species dataset
> are explicitly and individually described in this example, but it is
> sufficient in this case I think.
>
> \name{varespec}
> \alias{varechem}
> \alias{varespec}
> \docType{data}
> \title{Vegetation and environment in lichen pastures}
> \usage{
>data(varechem)
>data(varespec)
> }
> \description{
>   The \code{varespec} data frame has 24 rows and 44 columns.  Columns
>   are estimated cover values of 44 species.  The variable names are
>   formed from the scientific names, and are self explanatory for anybody
>   familiar with the vegetation type.
> The \code{varechem} data frame has 24 rows and 14 columns, giving the
> soil characteristics of the very same sites as in the \code{varespec}
> data frame. The chemical measurements have obvious names.
> \code{Baresoil} gives the estimated cover of bare soil, \code{Humpdepth}
> the thickness of the humus layer.
>
> }
> 
>
> HTH
>
> G
>
> > > -roger
> > >
> > > Arne Henningsen wrote:
> > > > Hi,
> > > >
> > > > since nobody answered to my first message, I try to explain my
> > > > problem more clearly and more general this time:
> > > >
> > > > I have a data set in my R package "micEcon", which has many variables
> > > > (82). Therefore, I would like to avoid to describe all variables in
> > > > the "\format" section of the documentation (.Rd file). However, doing
> > > > this lets "R CMD check" complain about "data codoc mismatches"
> > > > (details see below). Is there a way to avoid the description of all
> > > > variables without getting a complaint from "R CMD check"?
> > > >
> > > > Thanks,
> > > > Arne
> > > >
> > > >
> > > > --  Forwarded Message  --
> > > >
> > > > Subject: Documenting data sets with many variables
> > > > Date: Friday 05 August 2005 14:03
> > > > From: Arne Henningsen <[EMAIL PROTECTED]>
> > > > To: R-help@stat.math.ethz.ch
> > > >
> > > > Hi,
> > > >
> > > > I extended the data set "Blanciforti86" that is included in my R
> > > > package "micEcon". For instance, I added consumer prices, annual
> > > > consumption expenditures and expenditure shares of eleven aggregate
> > > > commodity groups. The corresponding variables in the data frame are
> > > > called "pAgg1", "pAgg2", ..., "pAgg11", "xAgg1", "xAgg2", ...,
> > > > "xAgg11", "wAgg1", "wAgg2", ..., "wAgg11". To avoid to describe all
> > > > 33 items in the "\format" section of the documentation (.Rd file) I
> > > > wrote something like
> > > >
> > > > \format{
> > > >This data frame contains the following columns:
> > > >\describe{
> > > >   [ . . . ]
> > > >   \item{xAggX}{Expenditure on the aggregate commodity group X
> > > >  (in Millions of US-Dollars).}
> > > >   \item{pAggX}{Price index for the aggregate commodity group X
> > > >  (1972 = 100).}
> > > >   \item{wAggX}{Expenditure share of the aggregate commodity group
> > > > X.} [ . . . ]
> > > >}
> > > > }
> > > >
> > > > and explained the 11 aggregate commodity groups only once in a
> > > > different section (1=food, 2=clothing, ... ). However, "R CMD check"
> > > > now complains about "data codoc mismatches", e.g.
> > > >   Code: [...] pAgg1pAgg2 pAgg3  [...]
> > > >   Docs: [...] pAggX [...]
> > > >
> > > > Is there a way to avoid the description of all 33 items without
> > > > getting a complaint from "R CMD check"?
> > > >
> > > > Thanks,
> > > > Arne
> > > >
> > > > ---

-- 
Arne Henningsen
Department of Agricultural Economics
University of Kiel
Olshausenstr. 40
D-24098 Kiel (Germany)
Tel: +49-431-880 4445
Fax: +49-431-880 1397
[EMAIL PROTECTED]
http://www.uni-kiel.de/agrarpol/ahenningsen/

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] optimal hedge variance ratio

[Krishna]

Hi All,

I am trying to finding out what could be an optimal hedge variance
ratio between spot and futures markets, between whose the degree of
correlation is highly varying.

For some reasons the hedge time period cannot extend for more than 1
month. So just to get an hint, I have calculated monthly correlation
coefficients which are highly varying. I am copying the frequency
distribution of monthly correlation coefficient values
(karl-pearsons') to indicating the degree of volatility.

Corre 
range   Frequency   
-0.70   0.00%
-0.52   3.08%
-0.33   4.62%
0   7   10.77%
0.3 7   10.77%
0.5 6   9.23%
0.7 16  24.62%
0.9 15  23.08%
1   9   13.85%

Can someone throw light on which model to use and how to approach for
desiging a hedge model (estimate hedge variance ratio) in such a
scenario. Help requested at the earliest.

thanks for the attention and best rgds

snvk

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] RODBC and sqlColumns

[Ben Stabler]

Ok, I understand that.  Then, how do I get the columns for a table that
is housed in a schema?  And, second, why does the following not work (or
at least partially work).  It creates the new table in the X schema but
then does not populate the table (and returns a sqlColumns() error.

sqlSave(db,x,"X.test",T,F)

Thanks.

Ben Stabler
Project Manager
PTV America, Inc.
1128 NE 2nd St, Suite 204
Corvallis, OR 97330
541-754-6836 x205
541-754-6837 fax
www.ptvamerica.com



-Original Message-
From: Prof Brian Ripley [mailto:[EMAIL PROTECTED] 
Sent: Wednesday, August 17, 2005 12:59 AM
To: Ben Stabler
Cc: r-help@stat.math.ethz.ch
Subject: Re: [R] RODBC and sqlColumns


AFAIK "." is not a valid part of an SQL table name. I think the help
files 
are perfectly clear as to what is supported:

  sqtable: character: a database table name accessible from the
   connected dsn.

Why do you think "X.test" is a `database table name'?

On Tue, 16 Aug 2005, Ben Stabler wrote:

> I have a Postgres database that I am connecting to with the Postgres 
> ODBC driver on Windows XP in R 2.1.0.  In the database is a database 
> with two schemas (public and X).  With RODBC (1.1-4) , I can connect 
> to the database and get the tables with sqlTables(db).  I can query 
> tables in the schema with sqlQuery("SELECT * FROM X.test").  However, 
> I can't get the columns in table X.test with sqlColumns(db,"X.test") 
> //it returns
>
> Error in sqlColumns(db, "X.test") : 'X.test': table not found on 
> channel
>
> If I do
>
> sqlColumns(db, "test") it returns
> [1] TABLE_QUALIFIER   TABLE_OWNER   TABLE_NAMECOLUMN_NAME
> DATA_TYPE
> [6] TYPE_NAME PRECISION LENGTHSCALE
> RADIX
> [11] NULLABLE  REMARKS   COLUMN_DEF
SQL_DATA_TYPE
> SQL_DATETIME_SUB
> [16] CHAR_OCTET_LENGTH ORDINAL_POSITION  IS_NULLABLE
DISPLAY_SIZE
> FIELD_TYPE
> <0 rows> (or 0-length row.names)
>
> But there is no test table defined anywhere else but the X schema.  If

> I do sqlSave(db,aDataFrame,"X.test",T,F), it says test already 
> defined. If I change the aDataFrame to be different than the fields 
> actually in the data, then R starts to create a new table but returns
>
> Error in sqlColumns(db, "X.test") : 'X.test': table not found on 
> channel
>
> It seems to be having problems with what is returned by the 
> columns.since
>
> Error in sqlSave(db, aDataFrame, "X.test", T, F) :
>[RODBC] ERROR: Could not SQLExecDirectS1000 7 ERROR:  relation 
> "test" already exists
>
> but if I change the input table to be differentthen R can create 
> the table, but fails to populate it.  I checked the db in PgAdmin and 
> the table is created by the sqlSave call.  All this stuff works if I 
> don't use a schema "schema.table".  So it appears there is something 
> wrong in some place dealing with understanding the columns for tables 
> in schemas.
>
> Any ideas?  Any help would be much appreciated.  Thank you.

-- 
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] Fitting mixture model

[Sean Davis]

I would like to fit a gaussian mixture model to a vector with about 50,000
points.  I have tried using Mclust to do so, but 50,000 points requires more
memory than I have (and I am running with 4Gb).  Any other suggestions for
how to do so?  Oh, I don't know the number of components, but the number
will likely be less than 5 or 6.

Thanks,
Sean

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] Fitting mixture model

[Christian Hennig]

I don't know if and how exactly it works for your data but package flexmix
should also allow to fit Gaussian mixtures.

Christian

On Wed, 17 Aug 2005, Sean Davis wrote:

> I would like to fit a gaussian mixture model to a vector with about 50,000
> points.  I have tried using Mclust to do so, but 50,000 points requires more
> memory than I have (and I am running with 4Gb).  Any other suggestions for
> how to do so?  Oh, I don't know the number of components, but the number
> will likely be less than 5 or 6.
>
> Thanks,
> Sean
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>

*** NEW ADDRESS! ***
Christian Hennig
University College London, Department of Statistical Science
Gower St., London WC1E 6BT, phone +44 207 679 1698
[EMAIL PROTECTED], www.homepages.ucl.ac.uk/~ucakche

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] Labels on axes with log scales with lattice

[Deepayan Sarkar]

On 8/14/05, Jamieson Cobleigh <[EMAIL PROTECTED]> wrote:
> I using lattice to make some plots and I want to make the y-axis on
> some of these plots use a log scale.  In the following plot:
> 
> x <- 1:10
> y <- 2^x
> xyplot(log10(y) ~ x)
> 
> I get tick marks on the y-axis at 0.5, 1.0, 1.5, 2.0, 2.5, and 3.0.  I
> would rather have just 3 tick marks at 1.0, 2.0, and 3.0 but labeled
> 10, 100, and 1000.
> 
> I know this can be done using the "at" and "labels" parameters to the
> "x" parameter to the "scales" parameter to the "xyplot" command.
> 
> xyplot(log10(y) ~ x, scales=list(y=list(at=c(1, 2, 3), labels=c(10,
> 100, 1000
> 
> My problem is that I am making multiple plots and cannot set the
> labels on each plot individually.  I need to automate the computation
> of the "at" and "labels" parameters.  I think the "axTicks" command
> can compute the information I need to set "at" and "labels" correctly,
> but I am having trouble determining how to set its parameters to make
> it compute the information I need.  Perhaps "pretty" might work to,
> but "axTicks" seems better designed for handling logarithmic axes.
> 
> Does anyone have any suggestions?

The `right' way to do this is 

xyplot(y ~ x, scales = list(y = list(log = 10)))

Unfortunately, the labeling for this doesn't use axTicks, it takes the
easy way out by using labels of the form "10^2", "10^3", etc. This is
partly due to laziness on my part, and also the fact that axTicks
doesn't support all the features 'scales' claims to.

My intended `solution' to this (currently vapourware) is to allow the
user to specify a function to calculate tick positions and labels. In
principle, this could be useful for other transformations, e.g. sqrt
for rootograms. I haven't thought through what the API for this would
be like, and I don't know when I will get around to actually
implementing it.

Deepayan

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] Dots in models formulae

[Laurent Valdes]

OK thanks everybody.

At least I should have try to use terms.formula().
But this is interesting.

2005/8/16, Prof Brian Ripley <[EMAIL PROTECTED]>:
> 
> On Tue, 16 Aug 2005, Gavin Simpson wrote:
> 
> > I still can't find this documented for a formula (I found
> > update.formula, but the meaning of "." is different there, slightly, as
> > you indicate) - but it must be as I didn't imagine seeing it - or maybe
> > I did...
> 
> It _is_ in `An Introduction to R', both for lm() and update.formula().
> Since this meaning does depend on using terms.formula(), it is on that
> functions' help page.
> 
> --
> Brian D. Ripley, [EMAIL PROTECTED]
> Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
> University of Oxford, Tel: +44 1865 272861 (self)
> 1 South Parks Road, +44 1865 272866 (PA)
> Oxford OX1 3TG, UK Fax: +44 1865 272595
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html
> 



-- 
--~~ Toulouse, Grenoble, Auch, Arcachon, Béziers, Paris, 
Saragosse, Lévignac Sur Save, habitats naturel du Valdo. ~~--


[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] Classifying values in vector

[mark salsburg]

I have a vector of size 217 called "A".

the values of A are not sorted and range from 0 to 1 (normalized)

I am having difficulty writing a program to create a new vector "B" where

if A's value is 0< A <=0.333 then B is 0
if A's value is 0.333< A <=0.666 then B is 1
if A's value is 0.666< A <=1 then B is 2

so if A is

0.22
0.999
0.444
0

B would be

0
2
1
0

thank you,

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] Classifying values in vector

[Sundar Dorai-Raj]

mark salsburg wrote:
> I have a vector of size 217 called "A".
> 
> the values of A are not sorted and range from 0 to 1 (normalized)
> 
> I am having difficulty writing a program to create a new vector "B" where
> 
> if A's value is 0< A <=0.333 then B is 0
> if A's value is 0.333< A <=0.666 then B is 1
> if A's value is 0.666< A <=1 then B is 2
> 
> so if A is
> 
> 0.22
> 0.999
> 0.444
> 0
> 
> B would be
> 
> 0
> 2
> 1
> 0
> 
> thank you,

Sounds like you are looking for ?cut:

 > A <- runif(217)
 > B <- cut(A, c(0, 1/3, 2/3, 1), labels = c(0, 1, 2))
 > # convert factor to numeric
 > as.numeric(levels(B)[B])

HTH,

--sundar

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] Classifying values in vector

[Sundar Dorai-Raj]

mark salsburg wrote:
> I have a vector of size 217 called "A".
> 
> the values of A are not sorted and range from 0 to 1 (normalized)
> 
> I am having difficulty writing a program to create a new vector "B" where
> 
> if A's value is 0< A <=0.333 then B is 0
> if A's value is 0.333< A <=0.666 then B is 1
> if A's value is 0.666< A <=1 then B is 2
> 
> so if A is
> 
> 0.22
> 0.999
> 0.444
> 0
> 
> B would be
> 
> 0
> 2
> 1
> 0
> 
> thank you,
> 

Didn't look closely at your example. If "A" can == zero then you need to 
include "include.lowest = TRUE" in your call to cut. This would mean 
your breaks are:

if A's value is 0 <= A <=0.333 then B is 0

which is different from what you have written above.


--sundar

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] Concurrence Matrix

[mark salsburg]

Thank you for all the help.

Can someone refer me to a function that can help with the creation of
a concurrence matrix?

I have two asymmetric vectors (C is length 217 and D is length 16063)
that contain values of 0,1,2

I want to create a matrix E (where rows are D and columns are C)  that
ouputs a score where

if C and D are the same number then E has value 1
if C and D are different numbers (i.e. 0 and 1, or 2 and 1) then E has a value 0

so for example:

so if C is c(0,0,2,1,0,0)
and  D is c(0,0,1,1,1,1,1,1)

E would look like:

  C  0  0  2  1  0  0

D0  1  1  0  0  0  0  
  0  1  1  0  0  1  1
  1  0  0  0  1  0  0
  1  0  0  0  1  0  0
  1  0  0  0  1  0  0 
  1  0  0  0  1  0  0
  1  0  0  0  1  0  0 
  1  0  0  0  1  0  0

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] accesing slots of S4 class in C code

[Aniko Szabo]

I am trying to use a custom S4 object in my C code and I cannot get the
access to its slots working.

 

The following is a toy example, but even this crashes. 

 

In R I have:

 

setClass("pd", representation(data="numeric"))

x <- new("pd", data=1:5)

 

test <- function(pd.obj) {

  res <- .C("TestFun", pd.obj)

  res}

 

test(x)

 

(Of couse I load the DLL as well.)

 

The corresponding C file has:

 

SEXP TestFun(SEXP pd_obj)

{

 SEXP t=R_NilValue;

 PROTECT(t = GET_SLOT(pd_obj, install("data")));

 UNPROTECT(1);

 return(t);

}

 

 

What I hope to get as a result is the (1:5) vector. 

In the long term, the vector will be a multi-dimensional array and I
will want to do calculations using its contents in the C program.

 

Thanks for any help,

 

Aniko


Huntsman Cancer Institute wishes to promote open communication while protecting 
confidential and/or privileged information.  If you have received this message 
in error, please inform the sender and delete all copies.


[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] Concurrence Matrix

[Berton Gunter]

?outer
Use "!=" for FUN and convert to numeric with as.numeric() if you like.

-- Bert Gunter
Genentech Non-Clinical Statistics
South San Francisco, CA
 
"The business of the statistician is to catalyze the scientific learning
process."  - George E. P. Box
 
 

> -Original Message-
> From: [EMAIL PROTECTED] 
> [mailto:[EMAIL PROTECTED] On Behalf Of mark salsburg
> Sent: Wednesday, August 17, 2005 10:26 AM
> To: R-help@stat.math.ethz.ch
> Subject: [R] Concurrence Matrix
> 
> Thank you for all the help.
> 
> Can someone refer me to a function that can help with the creation of
> a concurrence matrix?
> 
> I have two asymmetric vectors (C is length 217 and D is length 16063)
> that contain values of 0,1,2
> 
> I want to create a matrix E (where rows are D and columns are C)  that
> ouputs a score where
> 
> if C and D are the same number then E has value 1
> if C and D are different numbers (i.e. 0 and 1, or 2 and 1) 
> then E has a value 0
> 
> so for example:
> 
> so if C is c(0,0,2,1,0,0)
> and  D is c(0,0,1,1,1,1,1,1)
> 
> E would look like:
> 
>   C  0  0  2  1  0  0
> 
> D0  1  1  0  0  0  0  
>   0  1  1  0  0  1  1
>   1  0  0  0  1  0  0
>   1  0  0  0  1  0  0
>   1  0  0  0  1  0  0 
>   1  0  0  0  1  0  0
>   1  0  0  0  1  0  0 
>   1  0  0  0  1  0  0
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] Concurrence Matrix

[Ted Harding]

On 17-Aug-05 mark salsburg wrote:
> Thank you for all the help.
> 
> Can someone refer me to a function that can help with the creation of
> a concurrence matrix?
> 
> I have two asymmetric vectors (C is length 217 and D is length 16063)
> that contain values of 0,1,2
> 
> I want to create a matrix E (where rows are D and columns are C)  that
> ouputs a score where
> 
> if C and D are the same number then E has value 1
> if C and D are different numbers (i.e. 0 and 1, or 2 and 1) then E has
> a value 0
> 
> so for example:
> 
> so if C is c(0,0,2,1,0,0)
> and  D is c(0,0,1,1,1,1,1,1)
> 
> E would look like:
> 
>   C  0  0  2  1  0  0
> 
> D0  1  1  0  0  0  0  
>   0  1  1  0  0  1  1
>   1  0  0  0  1  0  0
>   1  0  0  0  1  0  0
>   1  0  0  0  1  0  0 
>   1  0  0  0  1  0  0
>   1  0  0  0  1  0  0 
>   1  0  0  0  1  0  0

It looks as though

  1*outer(D,C,"==")

does what you want (and shows up that E[1,5] and E[1,6] seem
to be wrong in your example, since D[1] = C[5] = 0 and
D[1] = C[6] = 0).

Best wishes,
Ted.



E-Mail: (Ted Harding) <[EMAIL PROTECTED]>
Fax-to-email: +44 (0)870 094 0861
Date: 17-Aug-05   Time: 18:56:49
-- XFMail --

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] RODBC and sqlColumns

[Andreas Hary]

Try append = F, that works for me.

A


- Original Message - 
From: "Ben Stabler" <[EMAIL PROTECTED]>
To: 
Sent: Wednesday, August 17, 2005 3:52 PM
Subject: Re: [R] RODBC and sqlColumns


> Ok, I understand that.  Then, how do I get the columns for a table that
> is housed in a schema?  And, second, why does the following not work (or
> at least partially work).  It creates the new table in the X schema but
> then does not populate the table (and returns a sqlColumns() error.
>
> sqlSave(db,x,"X.test",T,F)
>
> Thanks.
>
> Ben Stabler
> Project Manager
> PTV America, Inc.
> 1128 NE 2nd St, Suite 204
> Corvallis, OR 97330
> 541-754-6836 x205
> 541-754-6837 fax
> www.ptvamerica.com
>
>
>
> -Original Message-
> From: Prof Brian Ripley [mailto:[EMAIL PROTECTED]
> Sent: Wednesday, August 17, 2005 12:59 AM
> To: Ben Stabler
> Cc: r-help@stat.math.ethz.ch
> Subject: Re: [R] RODBC and sqlColumns
>
>
> AFAIK "." is not a valid part of an SQL table name. I think the help
> files
> are perfectly clear as to what is supported:
>
>  sqtable: character: a database table name accessible from the
>   connected dsn.
>
> Why do you think "X.test" is a `database table name'?
>
> On Tue, 16 Aug 2005, Ben Stabler wrote:
>
>> I have a Postgres database that I am connecting to with the Postgres
>> ODBC driver on Windows XP in R 2.1.0.  In the database is a database
>> with two schemas (public and X).  With RODBC (1.1-4) , I can connect
>> to the database and get the tables with sqlTables(db).  I can query
>> tables in the schema with sqlQuery("SELECT * FROM X.test").  However,
>> I can't get the columns in table X.test with sqlColumns(db,"X.test")
>> //it returns
>>
>> Error in sqlColumns(db, "X.test") : 'X.test': table not found on
>> channel
>>
>> If I do
>>
>> sqlColumns(db, "test") it returns
>> [1] TABLE_QUALIFIER   TABLE_OWNER   TABLE_NAMECOLUMN_NAME
>> DATA_TYPE
>> [6] TYPE_NAME PRECISION LENGTHSCALE
>> RADIX
>> [11] NULLABLE  REMARKS   COLUMN_DEF
> SQL_DATA_TYPE
>> SQL_DATETIME_SUB
>> [16] CHAR_OCTET_LENGTH ORDINAL_POSITION  IS_NULLABLE
> DISPLAY_SIZE
>> FIELD_TYPE
>> <0 rows> (or 0-length row.names)
>>
>> But there is no test table defined anywhere else but the X schema.  If
>
>> I do sqlSave(db,aDataFrame,"X.test",T,F), it says test already
>> defined. If I change the aDataFrame to be different than the fields
>> actually in the data, then R starts to create a new table but returns
>>
>> Error in sqlColumns(db, "X.test") : 'X.test': table not found on
>> channel
>>
>> It seems to be having problems with what is returned by the
>> columns.since
>>
>> Error in sqlSave(db, aDataFrame, "X.test", T, F) :
>>[RODBC] ERROR: Could not SQLExecDirectS1000 7 ERROR:  relation
>> "test" already exists
>>
>> but if I change the input table to be differentthen R can create
>> the table, but fails to populate it.  I checked the db in PgAdmin and
>> the table is created by the sqlSave call.  All this stuff works if I
>> don't use a schema "schema.table".  So it appears there is something
>> wrong in some place dealing with understanding the columns for tables
>> in schemas.
>>
>> Any ideas?  Any help would be much appreciated.  Thank you.
>
> -- 
> Brian D. Ripley,  [EMAIL PROTECTED]
> Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
> University of Oxford, Tel:  +44 1865 272861 (self)
> 1 South Parks Road, +44 1865 272866 (PA)
> Oxford OX1 3TG, UKFax:  +44 1865 272595
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] Raw data type transformations

[dhinds]

Not sure if this belongs here or on the devel list:

I've needed a more efficient way to manipulate raw binary data in R,
with more than the minimal raw transformation functions in the base
package.  So I've written a small package in C that effectively lets
me cast back and forth between raw vectors and other data types.

I've implemented four functions: rawToHex, hexToRaw, readRaw, and
writeRaw.  The first two convert between raw vectors and hex character
strings.  I use these to work around a limitation of ROracle not
directly supporting raw data types, but they may have other uses.
readRaw and writeRaw are similar to readBin and writeBin, but operate
on raw vectors in place of binary connections, and I took a few
shortcuts in my implementation (I handle a subset of data types and
didn't implement byte swapping).

[aside: Rinternals.h does not export RAW() for packages]

I first thought there might be a way to do this via connections but
didn't want to actually read and write files just to map between data
types.  One option would be to implement a rawConnection() analog of
the current textConnection(), which would associate raw vectors with
binary connections, and allow readBin()/writeBin().

Looking at the text connection code, it seems to me that it could be
extended to handle this fairly easily.  For input connections, it
appears that all that needs to be done is to add support for filling
the internal buffer from a raw vector (a few lines of code), and
adding a text_read() function (again just a few lines).  For output,
text_write() is similarly simple, and a few lines of code would need
to be added to text_vfprintf() to handle copying to a raw vector in
place of a character vector.

The text_* functions would then probably be better named internal_*
since they would handle internal binary as well as internal text
connections.

Does this sound like it would be a useful capability?

-- David Hinds

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] Copying rows from a matrix using a vector of indices

[Martin Lam]

Hi,

I am trying to use a vector of indices to select some
rows from a matrix. But before I can do that I somehow
need to convert 'combinations' into a list, since
'mode(combinations)' says it's 'numerical'. Any idea
how I can do that?

library("combinat")

combinations <- t(combn(8,2))

indices <- c(sample(1:length(combinations),10))

# convert
???

subset <- combinations[indices]

Thanks in advance,

Martin

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] trouble with wilcox.test

[Greg Hather]

I'm having trouble with the wilcox.test command in R. To demonstrate the 
anomalous behavior of wilcox.test, consider 

> wilcox.test(c(1.5,5.5), c(1:1), exact = F)$p.value
[1] 0.01438390
> wilcox.test(c(1.5,5.5), c(1:1), exact = T)$p.value
[1] 6.39808e-07 (this calculation takes noticeably longer). 
> wilcox.test(c(1.5,5.5), c(1:2), exact = T)$p.value
(R closes/crashes)


I believe that wilcox.test(c(1.5,5.5), c(1:1), exact = F)$p.value yields a 
bad result because of the normal approximation which R uses when exact = F. Any 
suggestions for how to compute wilcox.test(c(1.5,5.5), c(1:2), exact = 
T)$p.value? 

--Greg

[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] do glm with two data sets

[Hu, Ying (NIH/NCI)]

I have two data sets:
File1.txt: 
Name id1   id2   id3   ...
N10 1 0 ...
N20 1 1 ...
N31 1 -1...
...
 
File2.txt:
Group id1   id2   id3   ...
G1   1.22 1.34 2.44 ...
G2   2.33 2.56 2.56 ...
G3   1.56 1.99 1.46 ...
...
I like to do:
x1<-c(0,1,0,...)
y1<-c(1.22,1.34, 2.44, ...)
z1<-data.frame(x,y)
summary(glm(y1~x1,data=z1)
 
But I do the same thing by inputting the data sets from the two files
e <- read.table("file1.txt", header=TRUE,row.names=1)
g <- read.table("file2.txt", header=TRUE,row.names=1)
e1<-exp[1,]
g1<-geno[1,]
d1<-data.frame(g, e)
summary(glm(e1 ~ g1, data=d1))
 
the error message is 
Error in model.frame(formula, rownames, variables, varnames, extras,
extranames,  : 
invalid variable type
Execution halted
 
Thanks in advance,
 
Ying

[[alternative HTML version deleted]]

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] Copying rows from a matrix using a vector of indices

[Marc Schwartz (via MN)]

On Wed, 2005-08-17 at 12:02 -0700, Martin Lam wrote:
> Hi,
> 
> I am trying to use a vector of indices to select some
> rows from a matrix. But before I can do that I somehow
> need to convert 'combinations' into a list, since
> 'mode(combinations)' says it's 'numerical'. Any idea
> how I can do that?
> 
> library("combinat")
> 
> combinations <- t(combn(8,2))
> 
> indices <- c(sample(1:length(combinations),10))
> 
> # convert
> ???
> 
> subset <- combinations[indices]
> 
> Thanks in advance,
> 
> Martin


Your goal is not entirely clear.  Given your code above, you end up
with:

> library("combinat")
> 
> combinations <- t(combn(8,2))
> combinations
  [,1] [,2]
 [1,]12
 [2,]13
 [3,]14
 [4,]15
 [5,]16
 [6,]17
 [7,]18
 [8,]23
 [9,]24
[10,]25
[11,]26
[12,]27
[13,]28
[14,]34
[15,]35
[16,]36
[17,]37
[18,]38
[19,]45
[20,]46
[21,]47
[22,]48
[23,]56
[24,]57
[25,]58
[26,]67
[27,]68
[28,]78

combinations is a 28 x 2 matrix (or a 1d vector of length 56).

Your use of sample() with '1:length(combinations)' is the same as 1:56.

> length(combinations)
[1] 56

Hence, you end up with:

> set.seed(128)
> indices <- sample(1:length(combinations), 10)

> indices
 [1] 41 53 38 46 50 44 25 11 43  7


I suspect that you really want to use '1:nrow(combinations)', which
yields 1:28.

> nrow(combinations)
[1] 28


Thus,

> set.seed(128)
> indices <- sample(1:nrow(combinations), 10)

> indices
 [1] 21 26 18 22 23 20 11  5 27  3


Now, you can use:

> subset <- combinations[indices, ]
> subset
  [,1] [,2]
 [1,]47
 [2,]67
 [3,]38
 [4,]48
 [5,]56
 [6,]46
 [7,]26
 [8,]16
 [9,]68
[10,]14

which yields the rows from combinations, defined by 'indices'.

If correct, please review "An Introduction to R" and ?"[" for more
information on subsetting objects in R.

HTH,

Marc Schwartz

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] plotting issue with timestamps

[Dhiren DSouza]

I have a dataset with transactions and a timestamp at which they occoured 
during a day.  The time stamp is in the format /MM/DD hh:mm:ss.  I would 
like to plot a timeseries of the transactions to see if there is a 
particular time in the day when there is a spike in transactions.  Ofcourse 
the /MM/DD can be dropped since I am monitoring activity for the day and 
the actual date is unimportant.

Can anyone give me some direction on this.  I could possibly build a 
frequency table but not sure how to plot against the timestamp in the format 
above.  Any help would be appreciated.

-Dhiren

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] newbie- how to add a row to a matrix

[Simone Gabbriellini]

dear list,
if I have a matrix

s<-matrix(1:5, ncol=5)

how can I add another row with other data to that matrix?

thank you,
simone

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] newbie- how to add a row to a matrix

[Sundar Dorai-Raj]

?rbind

Simone Gabbriellini wrote:
> dear list,
> if I have a matrix
> 
> s<-matrix(1:5, ncol=5)
> 
> how can I add another row with other data to that matrix?
> 
> thank you,
> simone
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] GLM/GAM and unobserved heterogeneity

[Kyle G. Lundstedt]

Hello,
 I'm interested in correcting for and measuring unobserved  
heterogeneity ("missing variables") using R.  In particular, I'm  
searching for a simple way to measure the amount of unobserved  
heterogeneity remaining in a series of increasingly complex models  
(adding additional variables to each new model) on the same data.
 I have a static database of 400,000 or so individual mortgage  
loans, each of which is observed monthly from origination (t=0) until  
termination (a binary yes/no variable).  In my update database, there  
are up to 60 months of observed data for each loan in the static  
database, and an individual loan has an "average life" of roughly 36  
months.
 Each loan has static covariates observed at origination, such as  
original loan amount and credit score, as well as time-varying  
covariates (TVC) such as age, interest rates, and house prices.   
Because these TVC change each month, I've constructed a modeling  
database that merges the static database with the update database.
 The resulting "loan-month" modeling database has one observation  
for every loan-month, and the static covariates remain the same for  
all loan-months for a given loan.  Thus, the modeling database has  
roughly 14.4 million loan-month records.  A loan is considered  
"active" as long as it has not yet terminated or been censored; my  
interest is in predicting termination.
 This type of data is often referred to as "event history" or  
"discrete hazard" data.  The standard R package to apply to such data  
is "survival", with which I could estimate a Cox proportional hazard  
model using coxph.  The advantage of such an approach is that  
unobserved heterogeneity is easily addressed using the "frailty" term.
 The disadvantages, at least for my purposes, are two-fold.   
First, my audience is unfamiliar with hazard models.  Second, my  
monthly data has many "ties" (many terminations in the same month),  
so I've been told that coxph won't work well on a large dataset with  
many ties.
 On the other hand, because the data is measured discretely each  
month, many references suggest applying generalized linear models  
(GLM, "logit"-type models) or even generalized addivitive models  
(GAM, "logit"-type models that incorporate nonlinearity in individual  
covariates).  The advantage to this approach is that GLM and GAM are  
readily available in R, and my audience is very familiar with logit- 
type models.
 The disadvantage, however, is that I am totally unfamiliar with  
ways to correct for and measure unobserved heterogeneity using GLM/ 
GAM-type models.  I've been told that unobserved heterogeneity in the  
hazard framework is analogous to random effects in the GLM/GAM  
framework, but there seem to be a number of R packages that address  
this issue in different ways.
 So, I'd greatly appreciate suggestions on a simple way to  
incorporate unobserved heterogeneity into a GLM/GAM-type model.  I'm  
not much of a statistician, so simple examples are always helpful.   
I'm also happy to track down specific article/book references, if  
folks think those might be of help.

Many thanks,
Kyle
---
kyle  at  hotmail . com
(email altered in obvious ways)

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] plotting issue with timestamps

[Gabor Grothendieck]

On 8/17/05, Dhiren DSouza <[EMAIL PROTECTED]> wrote:
> I have a dataset with transactions and a timestamp at which they occoured
> during a day.  The time stamp is in the format /MM/DD hh:mm:ss.  I would
> like to plot a timeseries of the transactions to see if there is a
> particular time in the day when there is a spike in transactions.  Ofcourse
> the /MM/DD can be dropped since I am monitoring activity for the day and
> the actual date is unimportant.
> 
> Can anyone give me some direction on this.  I could possibly build a
> frequency table but not sure how to plot against the timestamp in the format
> above.  Any help would be appreciated.


Lets say your times are in a character vector tt.  (If the dates
are there too use substring to remove them.) Then convert them
to times class (of library chron) which will represent them
internally as a fraction of a day where the day starts at 0 and 
goes to 1 so that 0.5 is noon.  Since we want a more meaningful
axis do not use the automatic plot axis (xaxt = "n") and also
restrict the xaxis to the internval 0:1.  So that we can get
all the hours to fix reduce the size of the axis labels to 70%.
Then redraw it yourself with the hours labelled and add some 
light grey vertical lines if you like.

library(chron)
tt <- c("10:11:12", "09:10:11", "01:02:03")
tt.times <- times(tt)
plot(density(tt.times), xaxt = "n", xlim = 0:1, xlab = "Hour")
axis(1, 0:24/24, 0:24, cex.axis = 0.7)
abline(v = 0:24/24, col = "lightgrey")

See RNews 4/1 Help Desk article for more.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] do glm with two data sets

[Sundar Dorai-Raj]

Hu, Ying (NIH/NCI) wrote:
> I have two data sets:
> File1.txt: 
> Name id1   id2   id3   ...
> N10 1 0 ...
> N20 1 1 ...
> N31 1 -1...
> ...
>  
> File2.txt:
> Group id1   id2   id3   ...
> G1   1.22 1.34 2.44 ...
> G2   2.33 2.56 2.56 ...
> G3   1.56 1.99 1.46 ...
> ...
> I like to do:
> x1<-c(0,1,0,...)
> y1<-c(1.22,1.34, 2.44, ...)
> z1<-data.frame(x,y)
> summary(glm(y1~x1,data=z1)
>  
> But I do the same thing by inputting the data sets from the two files
> e <- read.table("file1.txt", header=TRUE,row.names=1)
> g <- read.table("file2.txt", header=TRUE,row.names=1)
> e1<-exp[1,]
> g1<-geno[1,]
> d1<-data.frame(g, e)
> summary(glm(e1 ~ g1, data=d1))
>  
> the error message is 
> Error in model.frame(formula, rownames, variables, varnames, extras,
> extranames,  : 
> invalid variable type
> Execution halted
>  
> Thanks in advance,
>  
> Ying

You have several inconsistencies in your example, so it will be 
difficult to figure out what you are trying to accomplish.

 > e <- read.table("file1.txt", header=TRUE,row.names=1)
 > g <- read.table("file2.txt", header=TRUE,row.names=1)
 > e1<-exp[1,]

What's "exp"? Also it's dangerous to use an R function as a variable 
name. Most of the time R can tell the difference, but in some cases it 
cannot.

 > g1<-geno[1,]

What's "geno"?

 > d1<-data.frame(g, e)

d1 is now e and g cbind'ed together?

 > summary(glm(e1 ~ g1, data=d1))

Are "e1" and "g1" elements of "d1"? From what you've told us, I don't 
know where the error is occurring. Also, if you are having errors, you 
can more easily isolate the problem by doing:

fit <- glm(e1 ~ g1, data = d1)
summary(fit)

This will at least tell you the problem is in your call to "glm" and not 
"summary.glm".

--sundar

P.S. Please (re-)read the POSTING GUIDE. Most of the time you will 
figure out problems such as these on your own during the process of 
creating a reproducible example.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] newbie- how to add a row to a matrix

[Martin Lam]

# row bind
a <- matrix(1:5)
a
a <- rbind(a, 6)
a

# column bind
b <- matrix(1:5)
b
b <- cbind(b, 6:12)
b
b <- cbind(b, 13)
b

Hope this helps,

Martin


--- Simone Gabbriellini <[EMAIL PROTECTED]> wrote:

> dear list,
> if I have a matrix
> 
> s<-matrix(1:5, ncol=5)
> 
> how can I add another row with other data to that
> matrix?
> 
> thank you,
> simone
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] do glm with two data sets

[Gavin Simpson]

On Wed, 2005-08-17 at 17:22 -0500, Sundar Dorai-Raj wrote:
> 
> Hu, Ying (NIH/NCI) wrote:
> > I have two data sets:
> > File1.txt: 
> > Name id1   id2   id3   ...
> > N10 1 0 ...
> > N20 1 1 ...
> > N31 1 -1...
> > ...
> >  
> > File2.txt:
> > Group id1   id2   id3   ...
> > G1   1.22 1.34 2.44 ...
> > G2   2.33 2.56 2.56 ...
> > G3   1.56 1.99 1.46 ...
> > ...
> > I like to do:
> > x1<-c(0,1,0,...)
> > y1<-c(1.22,1.34, 2.44, ...)
> > z1<-data.frame(x,y)
> > summary(glm(y1~x1,data=z1)
> >  
> > But I do the same thing by inputting the data sets from the two files
> > e <- read.table("file1.txt", header=TRUE,row.names=1)
> > g <- read.table("file2.txt", header=TRUE,row.names=1)
> > e1<-exp[1,]
> > g1<-geno[1,]
> > d1<-data.frame(g, e)
> > summary(glm(e1 ~ g1, data=d1))
> >  
> > the error message is 
> > Error in model.frame(formula, rownames, variables, varnames, extras,
> > extranames,  : 
> > invalid variable type
> > Execution halted
> >  
> > Thanks in advance,
> >  
> > Ying

Hi Ying,

That error message is likely caused by having a data.frame on the right
hand side (rhs) of the formula. You can't have a data.frame on the rhs
of a formula and g1 is still a data frame even if you only choose the
first row, e.g.:

dat <- as.data.frame(matrix(100, 10, 10))
class(dat[1, ])
[1] "data.frame"

You could try:

glm(e1 ~ ., data=g1[1, ])

and see if that works, but as Sundar notes, your post is a little
difficult to follow, so this may not do what you were trying to achieve.

HTH

Gav

> 
> You have several inconsistencies in your example, so it will be 
> difficult to figure out what you are trying to accomplish.
> 
>  > e <- read.table("file1.txt", header=TRUE,row.names=1)
>  > g <- read.table("file2.txt", header=TRUE,row.names=1)
>  > e1<-exp[1,]
> 
> What's "exp"? Also it's dangerous to use an R function as a variable 
> name. Most of the time R can tell the difference, but in some cases it 
> cannot.
> 
>  > g1<-geno[1,]
> 
> What's "geno"?
> 
>  > d1<-data.frame(g, e)
> 
> d1 is now e and g cbind'ed together?
> 
>  > summary(glm(e1 ~ g1, data=d1))
> 
> Are "e1" and "g1" elements of "d1"? From what you've told us, I don't 
> know where the error is occurring. Also, if you are having errors, you 
> can more easily isolate the problem by doing:
> 
> fit <- glm(e1 ~ g1, data = d1)
> summary(fit)
> 
> This will at least tell you the problem is in your call to "glm" and not 
> "summary.glm".
> 
> --sundar
> 
> P.S. Please (re-)read the POSTING GUIDE. Most of the time you will 
> figure out problems such as these on your own during the process of 
> creating a reproducible example.
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
Gavin Simpson [T] +44 (0)20 7679 5522
ENSIS Research Fellow [F] +44 (0)20 7679 7565
ENSIS Ltd. & ECRC [E] gavin.simpsonATNOSPAMucl.ac.uk
UCL Department of Geography   [W] http://www.ucl.ac.uk/~ucfagls/cv/
26 Bedford Way[W] http://www.ucl.ac.uk/~ucfagls/
London.  WC1H 0AP.
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] How to assess significance of random effect in lme4

[Spencer Graves]

  Is there some reason you are NOT using "anova", as in "Examples" 
section of "?lmer"?

  Permit me to summarize what I know about this, and I'll be pleased if 
someone else who thinks they know different would kindly enlighten me 
and others who might otherwise be misled if anything I say is 
inconsistent with the best literature available at the moment:

  1.  Doug Bates in his PhD dissertation and later in his book with Don 
Watts (1988) Nonlinear Regression Analysis and Its Applications (Wiley) 
split approximation errors in nonlinear least squares into "intrinsic 
curvature" and "parameter effects curvature".  He quantified these two 
problems in the context of roughly three dozen published examples, if my 
memory is correct, and found that in not quite all cases, the parameter 
effects were at least an order of magnitude greater than the intrinsic 
curvature.

  2.  In nonnormal situations, maximum likelihood is subject to more 
approximation error -- intrinsic curvature -- than "simple" nonlinear 
least squares.  However, I would expect this comparison to still be 
fairly accurate, even if the differences may not be quite as stark.

  3.  The traditional use of "standard errors" to judge statistical 
significance is subject to both intrinsic and parameter effects errors, 
while likelihood ratio procedures such as anova are subject only to the 
intrinsic curvature (assuming there are no substantive problems with 
nonconvergence).  Consequently, to judge statistical significance of an 
effect, anova is usually substantially better than the so-called Wald 
procedure using approximate standard errors, and is almost never worse. 
  If anyone knows of a case where this is NOT true, I'd like to know.

  4.  With parameters at a boundary as with variance components, the 
best procedure seems to double the p-value from a nested anova (unless 
the reported p-value is already large).  This is because the 
2*log(likelihood ratio) in such cases is roughly a 50-50 mixture of 0 
and chi-square(1) [if testing only 1 variance component parameter]. 
This is supported by a substantial amount of research, including 
simulations discussed in a chapter in Pinheiro and Bates (2000) 
Mixed-Effects Models in S and S-Plus (Springer).  The may be more 
accurate procedures available in the literature, but none so simple as 
this as far as I know.

  Comments?
  spencer graves
p.s.  It looks like [EMAIL PROTECTED] is a list containing vectors of length 29 
and 6 in your example.  I don't know what they are, but I don't see how 
they can be standard errors in the usual sense.

Doran, Harold wrote:

> These are the posterior variances of the random effects (I think more
> properly termed "empirical" posteriors).  Your model apparently includes
> three levels of random variation (commu, bcohort, residual). The first
> are the variances associated with your commu random effect and the
> second are the variances associated with the bcohort random effect.
> 
> Accessing either one would require
> 
> [EMAIL PROTECTED] or [EMAIL PROTECTED]
> 
> Obviously, replace "fm" with the name of your fitted model.
> 
> -Original Message-
> From: [EMAIL PROTECTED]
> [mailto:[EMAIL PROTECTED] On Behalf Of Shige Song
> Sent: Wednesday, August 17, 2005 7:50 AM
> To: r-help@stat.math.ethz.ch
> Subject: Re: [R] How to assess significance of random effect in lme4
> 
> Hi Harold,
> 
> Thanks for the reply. I looked at my outputs using str() as you
> suggested, here is the part you mentioned:
> 
>   ..@ bVar :List of 2
>   .. ..$ commu  : num [1, 1, 1:29] 5e-10 5e-10 5e-10 5e-10 5e-10 ...
>   .. ..$ bcohort: num [1, 1, 1:6] 1.05e-05 7.45e-06 6.53e-06 8.25e-06
> 7.11e-06 ...
> 
> where commu and bcohort are the two second-level units. Are these
> standard errors? Why the second vector contains a series of different
> numbers?
> 
> Thanks!
> 
> Shige
> 
> On 8/17/05, Doran, Harold <[EMAIL PROTECTED]> wrote:
> 
>> 
>>
>>You can extract the posterior variance of the random effect from the 
>>bVar slot of the fitted lmer model. It is not a hidden option, but a 
>>part of the fitted model. It just doesn't show up when you use
> 
> summary().
> 
>> 
>> Look at the structure of your object to see what is available using
> 
> str().
> 
>> 
>> However, your comment below seems to imply that it is incorrect for 
>>lmer to report SDs instead of the standard error, which is not true. 
>>That is a quantity of direct interest.
>> 
>> Other multilevel programs report the same exact statistics (for the 
>>most part). For instance, HLM reports the variances as well. If you 
>>want the posterior variance of an HLM model you need to extract it.
>>
>> 
>> 
>> -Original Message-
>> From:   [EMAIL PROTECTED] on behalf of
>>Shige Song
>> Sent:   Wed 8/17/2005 6:30 AM
>> To: r-help@stat.math.ethz.ch
>> Cc:
>> Subject:[R] How to assess significance of random effect in
> 
> lme4
> 

Re: [R] RMySQL not loading on Mac OS X

[Bill Northcott]

On 11/08/2005, at 8:00 PM, Georg Otto wrote:
> I have a problem loading RMySQL 0.5-5 on Mac OS 10.4.2 running R  
> 2.1.1.
>
> I installed RMySQL using:
>
> export PKG_CPPFLAGS="-I/usr/local/mysql/include"
> export PKG_LIBS="-L/usr/local/mysql/lib -lmysqlclient"
>
> R CMD INSTALL /Users/gwo/Desktop/RMySQL_0.5-5.tar.gz
>
>
> The installation seemed to work ok, but when I load RMySQL in R I get
> an error message:
>
>
>> library(RMySQL)
>>
> Error in dyn.load(x, as.logical(local), as.logical(now)) :
>  unable to load shared library '/Library/Frameworks/
> R.framework/Resources/library/RMySQL/libs/RMySQL.so':
>dlopen(/Library/Frameworks/R.framework/Resources/library/RMySQL/
> libs/RMySQL.so, 6): Symbol not found: _printf$LDBLStub
>Referenced from: /Library/Frameworks/R.framework/Resources/library/
> RMySQL/libs/RMySQL.so
>Expected in: flat namespace
> Error in library(RMySQL) : .First.lib failed for 'RMySQL'
>
> Any hint will be highly appreciated!
>
I notice this question never got a reply.

It would have been better asked on the sig-Mac list, but here are  
some pointers.

The
> libs/RMySQL.so, 6): Symbol not found: _printf$LDBLStub
is caused by not linking all the necessary system libraries.

Either
1.  there is an attempt to link objects and/or static libraries built  
with gcc-3.x/g77 with objects produced by gcc-4.x/gfortran.  This  
will not work.
or
2.  there is an attempt to link using ld or libtool rather than the  
gcc compiler driver which will ensure that appropriate system  
libraries are used.

FWIW I had no problem building it, but I was using an R package which  
I built from source.  So I know the same compiler was used throughout.

If you are using the R binary distribution, make sure you have run  
'sudo gcc_select 3.3' to get the right default compiler.

Bill Northcott

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] code a family of garch

[Nongluck Klibbua]

Dear R-helpers,

I was wondering if anyone has or knows someone who might have an
implementation
of algorithm for estimating garcht-t, egarch and gjr models. I try to
use Fseries but I don't know how to code these models.  
Thanks a million in advance,
Sincerely,
Nongluck

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] line/bar through median in lattice's "bwplot"?

[M. K.]

Is there a way to render a line through the median point in the boxplot
generated by the Lattice command "bwplot"?  The line basically bisects
the bar at the median point...

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] power of a matrix

[Spencer Graves]

Hi, Peter:

  I couldn't find "mexp" in the Matrix package, but I did find it in 
fMultivar and in Lindsey's rmutil.  These are different functions, but 
produced essentially the same answer for mexp(array(1:4, dim=c(2,2))). 
While hunting for that, I also also found reference by Doug Bates in a 
previous interchange on r-help to "a classic paper ... I would recommend 
reading":

  Moler C., van Loan C., (2003); _Nineteen dubious ways to compute
  the exponential of a matrix,  twenty-five years later_, SIAM
  Review 45, 3-49.

  This paper was cited in the help page for mexp in fMultivar but not 
in rmutil.

  spencer graves

Peter Dalgaard wrote:

> "Rau, Roland" <[EMAIL PROTECTED]> writes:
> 
> 
>>Thank you very much! Thanks also to the authors of this function,
>>Vincente Canto Cassola and Martin Maechler!
>>
>>This is exactly what I hoped for.
> 
> 
> 
>>>look at function ?mtx.exp() in the Malmig package, e.g.,
> 
> 
> Also, there is mexp() in the Matrix package. I'm not sure about the
> relative merits. mexp() is one of the less dubious implementations of
> matrix exponentials, but it does require to and from class "Matrix".
> mtx.exp is a bit unfortunately named as it appears to calculate matrix
> *powers* (which in this case is what you need).
> 

-- 
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

[EMAIL PROTECTED]
www.pdf.com 
Tel:  408-938-4420
Fax: 408-280-7915

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] power of a matrix

[Gabor Grothendieck]

Its expm.


On 8/17/05, Spencer Graves <[EMAIL PROTECTED]> wrote:
> Hi, Peter:
> 
>  I couldn't find "mexp" in the Matrix package, but I did find it in
> fMultivar and in Lindsey's rmutil.  These are different functions, but
> produced essentially the same answer for mexp(array(1:4, dim=c(2,2))).
> While hunting for that, I also also found reference by Doug Bates in a
> previous interchange on r-help to "a classic paper ... I would recommend
> reading":
> 
>  Moler C., van Loan C., (2003); _Nineteen dubious ways to compute
>  the exponential of a matrix,  twenty-five years later_, SIAM
>  Review 45, 3-49.
> 
>  This paper was cited in the help page for mexp in fMultivar but not
> in rmutil.
> 
>  spencer graves
> 
> Peter Dalgaard wrote:
> 
> > "Rau, Roland" <[EMAIL PROTECTED]> writes:
> >
> >
> >>Thank you very much! Thanks also to the authors of this function,
> >>Vincente Canto Cassola and Martin Maechler!
> >>
> >>This is exactly what I hoped for.
> >
> > 
> >
> >>>look at function ?mtx.exp() in the Malmig package, e.g.,
> >
> >
> > Also, there is mexp() in the Matrix package. I'm not sure about the
> > relative merits. mexp() is one of the less dubious implementations of
> > matrix exponentials, but it does require to and from class "Matrix".
> > mtx.exp is a bit unfortunately named as it appears to calculate matrix
> > *powers* (which in this case is what you need).
> >
> 
> --
> Spencer Graves, PhD
> Senior Development Engineer
> PDF Solutions, Inc.
> 333 West San Carlos Street Suite 700
> San Jose, CA 95110, USA
> 
> [EMAIL PROTECTED]
> www.pdf.com 
> Tel:  408-938-4420
> Fax: 408-280-7915
> 
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] Unexpected result of read.dbf

[Susumu Tanimura]

Hi there,

Unexpected result is given with read.dbf from foreign 0.8-9, which is
reproducible in R 2.1.1 with the following sample data:

 test.dbf 
KEYCODE,N,19,0
422010010
42201002101
42201002102
42201002103
42201002104
422010060
422010071
422010072
42201008001
42201008002
422010100
42201011001
42201011002
422010130
422010140
42201015001
42201015002
42201015003
422010180
422010190

In the example, the name of column is "KEYCODE" with numeric type, 19
digits, and no decimal.

The data was saved as CSV format for confirmation.

> dbf.test <- read.dbf("test.dbf")
> dbf.test[1:20,]
 [1] 422010010NANANANA 422010060 422010071
 [8] 422010072NANA 422010100NANA 422010130
[15] 422010140NANANA 422010180 422010190
> csv.test <- read.csv("test.csv")
> csv.test[1:20,]
 [1]   422010010 42201002101 42201002102 42201002103 42201002104   422010060
 [7]   422010071   422010072 42201008001 42201008002   422010100 42201011001
[13] 42201011002   422010130   422010140 42201015001 42201015002 42201015003
[19]   422010180   422010190

I wonder why I get NA from test.dbf.  I have this kind of trouble when
handling ESRI Shape files with maptools.  There is no choice to use
read.csv instead of read.dbf at executing read.shape.

--
Susumu Tanimura

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

[R] Problem with building R packages under Windows

[Dr L. Y Hin]

Dear all,

I am coming to the guru for advise here. I am a native user of Windows,
and S-plus, and the process of migrating my libraries from S to R has been
more than a painstaking task...

I am currently using R version 2.1.1 in Windows XP SP2.
I have read the "Writing R extensions", "the FAQ in R", and your valuable
document
"R for Windows Users", but still unable to compile a package in
R using Rcmd (likely stupidity on my part).

The followings are what I have done:

1. Downloaded and setup ActivePerl and the path (in control panel) of
c:\Perl\bin has been added automatically. I have downloaded version (5.8.7).

2. Downloaded and setup tools bundle from
www.stats.ox.ac.uk/pub/Rtools/tools.zip and
unzip them into c:\bin and added the path (in control panel) of c:\bin

3. Downloaded and setup the entire package of MikTex at
http://www.miktex.org/ and added
the path (in control panel) of C:\texmf\miktex\bin

4. Added the path (in control panel) referencing R (this version being
R2.0patched) which
is C:\Program Files\R\rw2000pat\bin

5. Since my codes are all in R scripts without compiled codes, I did not
install MinGW.

6. At DOS prompt at the preamble of the directory called foo (which
contained all subdirectories as created in R using
package.skeleton(name = "foo", list=cdgam.func)),
I typed the following:

C:\>Rcmd build --binary foo
Can't locate R/Dcf.pm in @INC <@INC contains: c:\share\perl c:/Perl/lib
c:/Perl/site/lib .> at C:/bin/build line 29.
BEGIN failed--compilation aborted at C:/bin/build line 29.

I'd be very grateful if anyone could kindly advise me on this matter.

Thanks.
Lin

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Re: [R] trouble with wilcox.test

[Prof Brian Ripley]

On Wed, 17 Aug 2005, Greg Hather wrote:

> I'm having trouble with the wilcox.test command in R.

Are you sure it is not the concepts that are giving 'trouble'?
What real problem are you trying to solve here?

> To demonstrate the anomalous behavior of wilcox.test, consider
>
>> wilcox.test(c(1.5,5.5), c(1:1), exact = F)$p.value
> [1] 0.01438390
>> wilcox.test(c(1.5,5.5), c(1:1), exact = T)$p.value
> [1] 6.39808e-07 (this calculation takes noticeably longer).
>> wilcox.test(c(1.5,5.5), c(1:2), exact = T)$p.value
> (R closes/crashes)
>
> I believe that wilcox.test(c(1.5,5.5), c(1:1), exact = F)$p.value 
> yields a bad result because of the normal approximation which R uses 
> when exact = F.

Expecting an approximation to be good in the tail for m=2 is pretty 
unrealistic.  But then so is believing the null hypothesis of a common 
*continuous* distribution.  Why worry about the distribution under a 
hypothesis that is patently false?

People often refer to this class of tests as `distribution-free', but they 
are not.  The Wilcoxon test is designed for power against shift 
alternatives, but here there appears to be a very large difference in 
spread.  So

> wilcox.test(5000+c(1.5,5.5), c(1:1), exact = T)$p.value
[1] 0.9989005

even though the two samples differ in important ways.


> Any suggestions for how to compute 
> wilcox.test(c(1.5,5.5), c(1:2), exact = T)$p.value?

I get (current R 2.1.1 on Linux)

> wilcox.test(c(1.5,5.5), c(1:2), exact = T)$p.value
[1] 1.59976e-07

and no crash.  So the suggestion is to use a machine adequate to the task, 
and that probably means an OS with adequate stack size.

>   [[alternative HTML version deleted]]

> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Please do heed it.  What version of R and what machine is this?  And do 
take note of the request about HTML mail.

-- 
Brian D. Ripley,  [EMAIL PROTECTED]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford, Tel:  +44 1865 272861 (self)
1 South Parks Road, +44 1865 272866 (PA)
Oxford OX1 3TG, UKFax:  +44 1865 272595

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html