Re: [R] I need help in seeing the code

[PIKAL Petr]

Hi

Anyway. AFAIK your code looks OK to me, provided you want find the number of 
rows without NA in each file.

Petr


> -Original Message-
> From: r-help-boun...@r-project.org [mailto:r-help-bounces@r-
> project.org] On Behalf Of Sarah Goslee
> Sent: Tuesday, July 22, 2014 8:22 AM
> To: Darius Mulia
> Cc: r-help@r-project.org
> Subject: Re: [R] I need help in seeing the code
>
> Hi Darius,
>
> This is the main R-help list. We don't know anything about your
> Coursera class, including what your code is supposed to do. Not only
> that, this list has a no homework policy.
>
> You need to use the discussion group associated with your course, or
> whatever materials it provides, to get assistance.
>
> Sarah
>
> On Tuesday, July 22, 2014, Darius Mulia  wrote:
>
> > Hi there,
> >
> > I am Darius and I am taking the R Programming course in Coursera. I
> > have a problem that I had spent so much looking for the problem. I
> > wrote my code and I believe that the code works perfectly fine
> because
> > it produces the result as what the course demanded. However, when I
> > tried to submit it, it says that my code is wrong. I do believe I
> make
> > mistake, but I cannot seem to find it. the code is as follow:
> >
> > Complete.R
> >
> > complete <- function(directory, id = 1:332) {
> >   file <- list.files(directory, full.names=TRUE)
> >   nobs <- c()
> >   for (i in id){
> > file1 <- read.csv(file[i])
> > nobs1 <- sum(complete.cases(file1))
> > nobs <- c(nobs, nobs1)
> > df <- data.frame(nobs)
> >   }
> >   return(data.frame(id,df))
> > }
> >
> >
> > Please give me a hint where I should look at.
> >
> > Thank you very much for your time and concern. I look forward hearing
> > back from you.
> >
> > Sincerely,
> >
> > Darius Mulia.
> >
> > [[alternative HTML version deleted]]
> >
> > __
> > R-help@r-project.org  mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> > http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
> >
>
>
> --
> Sarah Goslee
> http://www.stringpage.com
> http://www.sarahgoslee.com
> http://www.functionaldiversity.org
>
>   [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-
> guide.html
> and provide commented, minimal, self-contained, reproducible code.


Tento e-mail a jakékoliv k němu připojené dokumenty jsou důvěrné a jsou určeny 
pouze jeho adresátům.
Jestliže jste obdržel(a) tento e-mail omylem, informujte laskavě neprodleně 
jeho odesílatele. Obsah tohoto emailu i s přílohami a jeho kopie vymažte ze 
svého systému.
Nejste-li zamýšleným adresátem tohoto emailu, nejste oprávněni tento email 
jakkoliv užívat, rozšiřovat, kopírovat či zveřejňovat.
Odesílatel e-mailu neodpovídá za eventuální škodu způsobenou modifikacemi či 
zpožděním přenosu e-mailu.

V případě, že je tento e-mail součástí obchodního jednání:
- vyhrazuje si odesílatel právo ukončit kdykoliv jednání o uzavření smlouvy, a 
to z jakéhokoliv důvodu i bez uvedení důvodu.
- a obsahuje-li nabídku, je adresát oprávněn nabídku bezodkladně přijmout; 
Odesílatel tohoto e-mailu (nabídky) vylučuje přijetí nabídky ze strany příjemce 
s dodatkem či odchylkou.
- trvá odesílatel na tom, že příslušná smlouva je uzavřena teprve výslovným 
dosažením shody na všech jejích náležitostech.
- odesílatel tohoto emailu informuje, že není oprávněn uzavírat za společnost 
žádné smlouvy s výjimkou případů, kdy k tomu byl písemně zmocněn nebo písemně 
pověřen a takové pověření nebo plná moc byly adresátovi tohoto emailu případně 
osobě, kterou adresát zastupuje, předloženy nebo jejich existence je adresátovi 
či osobě jím zastoupené známá.

This e-mail and any documents attached to it may be confidential and are 
intended only for its intended recipients.
If you received this e-mail by mistake, please immediately inform its sender. 
Delete the contents of this e-mail with all attachments and its copies from 
your system.
If you are not the intended recipient of this e-mail, you are not authorized to 
use, disseminate, copy or disclose this e-mail in any manner.
The sender of this e-mail shall not be liable for any possible damage caused by 
modifications of the e-mail or by delay with transfer of the email.

In case that this e-mail forms part of business dealings:
- the sender reserves the right to end negotiations about entering into a 
contract in any time, for any reason, and without stating any reasoning.
- if the e-mail contains an offer, the recipient is entitled to immediately 
accept such offer; The sender of this e-mail (offer) excludes any acceptance of 
the offer on the part of the recipient containing any amendment or variation.
- the sende

Re: [R] I need help in seeing the code

[Sarah Goslee]

Hi Darius,

This is the main R-help list. We don't know anything about your Coursera
class, including what your code is supposed to do. Not only that, this list
has a no homework policy.

You need to use the discussion group associated with your course, or
whatever materials it provides, to get assistance.

Sarah

On Tuesday, July 22, 2014, Darius Mulia  wrote:

> Hi there,
>
> I am Darius and I am taking the R Programming course in Coursera. I have a
> problem that I had spent so much looking for the problem. I wrote my code
> and I believe that the code works perfectly fine because it produces the
> result as what the course demanded. However, when I tried to submit it, it
> says that my code is wrong. I do believe I make mistake, but I cannot seem
> to find it. the code is as follow:
>
> Complete.R
>
> complete <- function(directory, id = 1:332) {
>   file <- list.files(directory, full.names=TRUE)
>   nobs <- c()
>   for (i in id){
> file1 <- read.csv(file[i])
> nobs1 <- sum(complete.cases(file1))
> nobs <- c(nobs, nobs1)
> df <- data.frame(nobs)
>   }
>   return(data.frame(id,df))
> }
>
>
> Please give me a hint where I should look at.
>
> Thank you very much for your time and concern. I look forward hearing back
> from you.
>
> Sincerely,
>
> Darius Mulia.
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org  mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>


-- 
Sarah Goslee
http://www.stringpage.com
http://www.sarahgoslee.com
http://www.functionaldiversity.org

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] I need help in seeing the code

[Darius Mulia]

Hi there,

I am Darius and I am taking the R Programming course in Coursera. I have a
problem that I had spent so much looking for the problem. I wrote my code
and I believe that the code works perfectly fine because it produces the
result as what the course demanded. However, when I tried to submit it, it
says that my code is wrong. I do believe I make mistake, but I cannot seem
to find it. the code is as follow:

Complete.R

complete <- function(directory, id = 1:332) {
  file <- list.files(directory, full.names=TRUE)
  nobs <- c()
  for (i in id){
file1 <- read.csv(file[i])
nobs1 <- sum(complete.cases(file1))
nobs <- c(nobs, nobs1)
df <- data.frame(nobs)
  }
  return(data.frame(id,df))
}


Please give me a hint where I should look at.

Thank you very much for your time and concern. I look forward hearing back
from you.

Sincerely,

Darius Mulia.

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] 1st el of a list of vectors

[Hervé Pagès]

Hi Carol,

On 07/21/2014 09:10 PM, Richard M. Heiberger wrote:

l = list(c(1,2), c(3,5,6), c(7))

sapply(l, `[`, 1)


Using sapply() works but won't be very efficient if you have a very long
list. If you worry about efficiency, you can do the following (using the
IRanges package from Bioconductor):

  > library(IRanges)
  > eltlens <- elementLengths(l)
  > unlist(l, use.names=FALSE)[cumsum(eltlens) - eltlens + 1L]
  [1] 1 3 7

Only worth if the length of your list is > 10 though...

Cheers,
H.

PS: See http://bioconductor.org/packages/release/bioc/html/IRanges.html
for how to install the IRanges package.



On Mon, Jul 21, 2014 at 3:55 PM, carol white  wrote:

Hi,
If we have a list of vectors of different lengths, how is it possible to 
retrieve the first element of the vectors of the list?


l = list(c(1,2), c(3,5,6), c(7))

1,3,7 should be retrieved

Thanks

Carol
 [[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



--
Hervé Pagès

Program in Computational Biology
Division of Public Health Sciences
Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N, M1-B514
P.O. Box 19024
Seattle, WA 98109-1024

E-mail: hpa...@fhcrc.org
Phone:  (206) 667-5791
Fax:(206) 667-1319

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] odd, even indices of a vector

[Hervé Pagès]

Hi Carol,

On 07/21/2014 01:33 PM, carol white wrote:

Might be a trivial question but how to identify the odd and even indices of a 
vector?

x = c(1,z,w,2,6,7)

el of odd indices= 1,w,6
el of even indices= z,2,7


The easiest way is to subset your vector with c(TRUE, FALSE) to keep
only the odd indices:

  > letters[c(TRUE, FALSE)]
  [1] "a" "c" "e" "g" "i" "k" "m" "o" "q" "s" "u" "w" "y"

and with c(FALSE, TRUE) to keep only the even indices:

  > letters[c(FALSE, TRUE)]
  [1] "b" "d" "f" "h" "j" "l" "n" "p" "r" "t" "v" "x" "z"

Note that the above doesn't work properly with vector of length < 2.

A method that works independently of the length of the vector (using
the is.odd and is.even functions from the "identifying odd or even
number" post that you referred to):

  x[is.odd(seq_along(x))]

  x[is.even(seq_along(x))]

Hope this helps,
H.



given the def of odd and even in 
https://stat.ethz.ch/pipermail/r-help/2010-July/244299.html

should a loop be used?

for (i in 1: length(x))
if (is.odd(i)) print (i)


Carol

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



--
Hervé Pagès

Program in Computational Biology
Division of Public Health Sciences
Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N, M1-B514
P.O. Box 19024
Seattle, WA 98109-1024

E-mail: hpa...@fhcrc.org
Phone:  (206) 667-5791
Fax:(206) 667-1319

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] odd, even indices of a vector

[Bert Gunter]

Ms. White:

Unless I have seriously misjudged, you really really really need to go
through an R tutorial -- An Intro to R ships with R, but there are
many on the web -- before posting here further. You do not appear to
have made much of an effort to learn even the basics, and I consider
it unfair to post questions here on basic matters that a tutorial
would tell you about.

And, of course, my apologies if I have gotten it wrong.

Cheers,
Bert

Bert Gunter
Genentech Nonclinical Biostatistics
(650) 467-7374

"Data is not information. Information is not knowledge. And knowledge
is certainly not wisdom."
Clifford Stoll




On Mon, Jul 21, 2014 at 1:33 PM, carol white  wrote:
> Might be a trivial question but how to identify the odd and even indices of a 
> vector?
>
> x = c(1,z,w,2,6,7)
>
> el of odd indices= 1,w,6
> el of even indices= z,2,7
>
> given the def of odd and even in 
> https://stat.ethz.ch/pipermail/r-help/2010-July/244299.html
>
> should a loop be used?
>
> for (i in 1: length(x))
> if (is.odd(i)) print (i)
>
>
> Carol
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] 1st el of a list of vectors

[Richard M. Heiberger]

l = list(c(1,2), c(3,5,6), c(7))

sapply(l, `[`, 1)

On Mon, Jul 21, 2014 at 3:55 PM, carol white  wrote:
> Hi,
> If we have a list of vectors of different lengths, how is it possible to 
> retrieve the first element of the vectors of the list?
>
>
> l = list(c(1,2), c(3,5,6), c(7))
>
> 1,3,7 should be retrieved
>
> Thanks
>
> Carol
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] odd, even indices of a vector

[Richard M. Heiberger]

## these functions assume the argument is integer

odd <- function(x) x%%2 != 0

even <- function(x) x%%2 == 0

evenb <- function(x) !odd(x)

odd(1:10)
even(1:10)
evenb(1:10)

On Mon, Jul 21, 2014 at 4:33 PM, carol white  wrote:
> Might be a trivial question but how to identify the odd and even indices of a 
> vector?
>
> x = c(1,z,w,2,6,7)
>
> el of odd indices= 1,w,6
> el of even indices= z,2,7
>
> given the def of odd and even in 
> https://stat.ethz.ch/pipermail/r-help/2010-July/244299.html
>
> should a loop be used?
>
> for (i in 1: length(x))
> if (is.odd(i)) print (i)
>
>
> Carol
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Maximum likelihood estimation (stats4::mle)

[David Winsemius]

On Jul 21, 2014, at 12:10 PM, Ronald Kölpin wrote:

> Dear R-Community,
> 
> I'm trying to estimate the parameters of a probability distribution
> function by maximum likelihood estimation (using the stats4 function
> mle()) but can't seem to get it working.
> 
> For each unit of observation I have a pair of observations (a, r)
> which I assume (both) to be log-normal distributed (iid). Taking the
> log of both values I have (iid) normally distributed random variables
> and the likelihood function to be estimated is:
> 
> L = Product(F(x_i) - F(y_i), i=1..n)
> 
> where F is the Normal PDF and (x,y) := (log(a), log(r)). Taking the
> log and multiplying by -1 gives the negative loglikelihood

I don't see the need to multiply by -1. The log of any probability is of 
necessity less than (or possibly equal to) 0 since probabilities are bounded 
above by 1. So sums of these number will be negative which then allows you to 
maximize their sums.

> 
> l = Sum(log( F(x_i) - F(y_i) ), i=1..n)
> 
> However estimation by mle() produces the error "vmmin is not finite"


As I would have predicted. If one maximizes numbers that get larger as 
probabilities get small this is what would be expected.

-- 
David.

> and "NaN have been created" - even though put bound on the parameters
> mu and sigma (see code below).
> 
> 
> library("stats4")
> 
> gaps <- matrix(nrow=10, ncol=4, dimnames=list(c(1:10),c("r_i", "a_i",
> "log(r_i)", "log(a_i)")))
> gaps[,1] <- c(2.6, 1.4, 2.2, 2.9, 2.9, 1.7, 1.3, 1.7, 3.8, 4.5)
> gaps[,2] <- c(9.8, 20.5, 8.7, 7.2, 10.3, 11, 4.5, 5.2, 6.7, 7.6)
> gaps[,3] <- log(gaps[,1])
> gaps[,4] <- log(gaps[,2])
> 
> nll <- function(mu, sigma)
> {
>if(sigma >= 0 && mu >= 0)
>{
>-sum(log(pnorm(gaps[,3], mean=mu, sd=sigma) - pnorm(gaps[,4],
> mean=mu, sd=sigma)))
>}
>else
>{
>NA
>}
> }
> 
> fit <- mle(nll, start=list(mu=0, sigma=1), nobs=10)
> print(fit)
> 
> 
> To be honest, I'm stumped and don't quite know what the problem is...
> 
> Regards and Thanks,
> 
> Ronald Kölpin
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

David Winsemius
Alameda, CA, USA

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Help with SEM package - model significance

[Bernardo Santos]

Hi John,

Thanks for your reply (1 month later lol).
In fact maybe the point is that I do not understand exactly the role of latent 
variables (what they are, and how to define them in R) in SEM.
Do you have any suggestion of easy basic literature on SEM that can help me 
with that?
Most things I have read are old (and use some different statistics programs) 
and offer examples too far from ecology, than I had some difficulties to 
understand the method in general.

But, as far as I understood, SEM is like a simple multiple regression (linear 
model), but that takes into account the relation of different variables 
simultaneously, isn't that?

Thank you very much.
Best regards,

Bernardo



Em Segunda-feira, 16 de Junho de 2014 8:40, John Fox  
escreveu:
 


Dear Bernardo,

The df for the LR chisquare over-identification test come not from the number 
of observations, but from the difference between the number of observable 
variances and covariances, on the one hand, and free parameters to estimate, on 
the other. In your case, these numbers are equal, and so df = 0. The LR 
chisquare for a just-identified model is also necessarily 0: the model 
perfectly reproduces the covariational structure of the observed variables. 

R (and most statistical software) by default writes very small and very large 
numbers in scientific format. In your case, -2.873188e-13 = -2.87*10^-13, that 
is, 0 within rounding error. You can change the way numbers are printed with 
the R scipen option.

Some other observations:

(1) Your model is recursive and has no latent variables; you would get the same 
estimates from OLS regression using lm().

(2) For quite some time now, the sem package has included specifyEquations() as 
a more convenient way of specifying a model, in preference to specifyModel(). 
See ?specifyEquations.

(3) You don't have to specify the error variances directly; specifyEquations(), 
or specifyModel(), will supply them.

I hope this helps,
John



John Fox, Professor
McMaster University
Hamilton, Ontario, Canada
http://socserv.mcmaster.ca/jfox/
    
    
On Sun, 15 Jun 2014 20:15:31 -0700 (PDT)

> Dear all, 
> 
> I used "sem" function from the package SEM to fit a model. However, I cannot 
> say if the model is correspondent to the data or not (chisquare test).
> I used the commands:
> 
> model1 <- specifyModel()
> estadio -> compflora, a1, NA
> estadio -> compfauna, a2, NA
> estadio -> interacoesobs, a3, NA
> compflora -> compfauna, b1, NA
> compflora -> interacoesobs, b2, NA
> compfauna -> interacoesobs, c1, NA
> estadio <-> estadio, e1, NA
> compflora <-> compflora, e2, NA
> compfauna <-> compfauna, e3, NA
> interacoesobs <-> interacoesobs, e4, NA
> 
> sem1 <- sem(model1, cov.matrix, length(samples))
> summary(sem1)
> 
> and I got the result:
> 
> Model Chisquare =  -2.873188e-13   Df =  0 Pr(>Chisq) = NA AIC =  20 BIC =  
> -2.873188e-13 Normalized Residuals Min.   1st Qu.    Median      Mean   3rd 
> Qu.      Max. 
> 0.000e+00 0.000e+00 2.957e-16 3.193e-16 5.044e-16 8.141e-16  R-square for 
> Endogenous Variables compflora     compfauna interacoesobs  0.0657        
> 0.1056        0.2319  Parameter Estimates Estimate     Std Error    z value   
>  Pr(>|z|)                                    
> a1 3.027344e-01 1.665395e-01 1.81779316 6.909575e-02 compflora <--- estadio   
>        
> a2 2.189427e-01 1.767404e-01 1.23878105 2.154266e-01 compfauna <--- estadio   
>        
> a3 7.314192e-03 1.063613e-01 0.06876742 9.451748e-01 interacoesobs <--- 
> estadio      
> b1 2.422906e-01 1.496290e-01 1.61927587 1.053879e-01 compfauna <--- compflora 
>        
> b2 3.029933e-01 9.104901e-02 3.32780446 8.753328e-04 interacoesobs <--- 
> compflora    
> c1 4.863368e-02 8.638177e-02 0.56300857 5.734290e-01 interacoesobs <--- 
> compfauna    
> e1 6.918133e+04 1.427102e+04 4.84767986 1.249138e-06 estadio <--> estadio     
>        
> e2 9.018230e+04 1.860319e+04 4.84767986 1.249138e-06 compflora <--> compflora 
>        
> e3 9.489661e+04 1.957568e+04 4.84767986 1.249138e-06 compfauna <--> compfauna 
>        
> e4 3.328072e+04 6.865289e+03 4.84767986 1.249138e-06 interacoesobs <--> 
> interacoesobs Iterations =  0 
> 
> I understand the results, but I do not know how to interpret the first line 
> that tells me about the model:
> Model Chisquare =  -2.873188e-13   Df =  0 Pr(>Chisq) = NA
> 
> How can DF be zero, if the number of observations I used in sem funcition was 
> 48 and I have only 4 variables? What is the p value?
> 
> Thanks in advance.
> Bernardo Niebuhr
>     [[alternative HTML version deleted]]
> 
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] odd, even indices of a vector

[carol white]

Might be a trivial question but how to identify the odd and even indices of a 
vector?

x = c(1,z,w,2,6,7)

el of odd indices= 1,w,6
el of even indices= z,2,7

given the def of odd and even in 
https://stat.ethz.ch/pipermail/r-help/2010-July/244299.html

should a loop be used?

for (i in 1: length(x))
if (is.odd(i)) print (i)


Carol

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Maximum likelihood estimation (stats4::mle)

[Ronald Kölpin]

Dear R-Community,

I'm trying to estimate the parameters of a probability distribution
function by maximum likelihood estimation (using the stats4 function
mle()) but can't seem to get it working.

For each unit of observation I have a pair of observations (a, r)
which I assume (both) to be log-normal distributed (iid). Taking the
log of both values I have (iid) normally distributed random variables
and the likelihood function to be estimated is:

L = Product(F(x_i) - F(y_i), i=1..n)

where F is the Normal PDF and (x,y) := (log(a), log(r)). Taking the
log and multiplying by -1 gives the negative loglikelihood

l = Sum(log( F(x_i) - F(y_i) ), i=1..n)

However estimation by mle() produces the error "vmmin is not finite"
and "NaN have been created" - even though put bound on the parameters
mu and sigma (see code below).


library("stats4")

gaps <- matrix(nrow=10, ncol=4, dimnames=list(c(1:10),c("r_i", "a_i",
"log(r_i)", "log(a_i)")))
gaps[,1] <- c(2.6, 1.4, 2.2, 2.9, 2.9, 1.7, 1.3, 1.7, 3.8, 4.5)
gaps[,2] <- c(9.8, 20.5, 8.7, 7.2, 10.3, 11, 4.5, 5.2, 6.7, 7.6)
gaps[,3] <- log(gaps[,1])
gaps[,4] <- log(gaps[,2])

nll <- function(mu, sigma)
{
if(sigma >= 0 && mu >= 0)
{
-sum(log(pnorm(gaps[,3], mean=mu, sd=sigma) - pnorm(gaps[,4],
mean=mu, sd=sigma)))
}
else
{
NA
}
}

fit <- mle(nll, start=list(mu=0, sigma=1), nobs=10)
print(fit)


To be honest, I'm stumped and don't quite know what the problem is...

Regards and Thanks,

Ronald Kölpin

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] standard error of survfit.coxph()

[array chip]

Dear Terry/All,

I was trying to use your explanation of the standard error estimate from 
survfit.coxph() to verify the standard error estimates for the method of 
log(log(S)), but couldn't get the estimates correct. Here is an example using 
the lung dataset:

> fit<-coxph(Surv(time,status)~wt.loss,lung)
> surv<-survfit(fit,newdata=lung[2:5,])$surv

> logstd<-survfit(fit,newdata=lung[2:5,])$std.err  
> logstd.time<-survfit(fit,newdata=lung[2:5,])$time


    ## std error of cumulative hazard at time=60

> logstd<-logstd[logstd.time==60,]
> logstd
[1] 0.01965858 0.01965858 0.01941871 0.01969248

    ## survival S at months 60

> surv<-surv[logstd.time==60,]  
> surv
    2 3 4 5 
0.9249238 0.9249238 0.9253038 0.9263392

Since survfit()$std.err was for cumulative hazard H=log(S), thus based on your 
explanation below, the standard error for log-log method is: (1/S^2)var(H)

> loglogstd<-sqrt(1/surv^2*(logstd^2))
> loglogstd

 2  3  4  5 
0.02125427 0.02125427 0.02098631 0.02125839

    ## the upper confidence interval should be
> exp(-exp(log(-log(surv))-qnorm(0.975)*loglogstd))
    2
 3 4 5 
0.9278737 0.9278737 0.9282031 0.9292363

But this is different from the output using summary.survfit:

> summary(survfit(fit,newdata=lung[2:5,],conf.type='log-log'),times=60)$upper[1,]
    2 3 4 5 
0.9534814 0.9534814 0.9535646 0.9548478

Can you please suggest what I did wrong in the calculation?

Thanks very much!

John




 From: "Therneau, Terry M., Ph.D." 

Sent: Monday, June 30, 2014 6:04 AM
Subject: Re:  standard error of survfit.coxph()


1. The computations "behind the scenes" produce the variance of the cumulative 
hazard. 
This is true for both an ordinary Kaplan-Meier and a Cox model.  
Transformations to other 
scales are done using simple Taylor series.

   H = cumulative hazard = log(S);  S=survival
   var(H) = var(log(S))  = the starting point
   S = exp(log(S)), so  var(S) is approx [deriv of exp(x)]^2 * var(log(S)) = 
S^2 var(H)
   var(log(log(S)) is approx (1/S^2) var(H)

2. At the time it was written, summary.survfit was used only for printing out 
the survival 
curve at selected times, and the audience for the printout wanted std(S).   
True, that was 
20 years ago, but I don't recall anyone ever asking for summary to do anything 
else.  Your 
request is not a bad idea.
   Note however that the primary impact of using log(S) or S or log(log(S)) 
scale is is on 
the confidence intervals, and they do appear per request in the summary output.

Terry T.


On 06/28/2014 05:00 AM, r-help-requ...@r-project.org wrote:
> Message: 9
> Date: Fri, 27 Jun 2014 12:39:29 -0700

> To:"r-help@r-project.org"  
> Subject: [R] standard error of survfit.coxph()
> Message-ID:
>     <1403897969.91269.yahoomail...@web122906.mail.ne1.yahoo.com>
> Content-Type: text/plain
>
> Hi, can anyone help me to understand the standard errors printed in the 
> output of survfit.coxph()?
>
> time<-sample(1:15,100,replace=T)
>
> status<-as.numeric(runif(100,0,1)<0.2)
> x<-rnorm(100,10,2)
>
> fit<-coxph(Surv(time,status)~x)
> ??? ### method 1
>
> survfit(fit, newdata=data.frame(time=time,status=status,x=x)[1:5,], 
> conf.type='log')$std.err
>
>  [,1]??? [,2]??? [,3]??? [,4]?? [,5]
> ?[1,] 0.0 0.0 0.0 0.0 0.
> ?[2,] 0.008627644 0.008567253 0.008773699 0.009354788 0.01481819
> ?[3,] 0.008627644 0.008567253 0.008773699 0.009354788 0.01481819
> ?[4,] 0.013800603 0.013767977 0.013889971 0.014379928 0.02353371
> ?[5,] 0.013800603 0.013767977 0.013889971 0.014379928 0.02353371
> ?[6,] 0.013800603 0.013767977 0.013889971 0.014379928 0.02353371
> ?[7,] 0.030226811 0.030423883 0.029806263 0.028918817 0.05191161
> ?[8,] 0.030226811 0.030423883 0.029806263 0.028918817 0.05191161
> ?[9,] 0.036852571 0.037159980 0.036186931 0.034645002 0.06485394
> [10,] 0.044181716 0.044621159 0.043221145 0.040872939 0.07931028
> [11,] 0.044181716 0.044621159 0.043221145 0.040872939 0.07931028
> [12,] 0.055452631 0.056018832 0.054236881 0.051586391 0.10800413
> [13,] 0.070665160 0.071363749 0.069208056 0.066655730 0.14976433
> [14,] 0.124140400 0.125564637 0.121281571 0.118002021 0.30971860
> [15,] 0.173132357 0.175309455 0.168821266 0.164860523 0.46393111
>
> survfit(fit, newdata=data.frame(time=time,status=status,x=x)[1:5,], 
> conf.type='log')$time
> ?[1]? 1? 2? 3? 4? 5? 6? 7? 8? 9 10 11 12 13 14 15
>
> ??? ### method 2
>
> summary(survfit(fit, newdata=data.frame(time=time,status=status,x=x)[1:5,], 
> conf.type='log'),time=10)$std.err
>
> ? 1? 2? 3? 4? 5
> [1,] 0.04061384 0.04106186 0.03963184 0.03715246 0.06867532
>
> By reading the help of ?survfit.object and ?summary.survfit, the standard 
> error provided in the output of method 1 (survfit()) was for cumulative 
> hazard-log(survival), while the standard

[R] 1st el of a list of vectors

[carol white]

Hi,
If we have a list of vectors of different lengths, how is it possible to 
retrieve the first element of the vectors of the list? 


l = list(c(1,2), c(3,5,6), c(7))

1,3,7 should be retrieved

Thanks

Carol
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] duplicated rows of a matrix

[carol white]

I need that duplicated indicate all row indices (occurences) that are 
duplicated. In your example, rows 2,3,4,5,8

Thanks.

Carol



On Monday, July 21, 2014 9:17 PM, William Dunlap  wrote:
 


Can you give an example of duplicated() not working on the rows of a matrix?

Here is an example where it does work:
> m <- cbind(c(a=1,b=2,c=3,d=2,e=3,f=4,g=1,h=1), c(11,13,11,13,11,13,13,11))
> class(m)
[1] "matrix"
> m
  [,1] [,2]
a    1   11
b    2   13
c    3   11
d    2   13
e    3   11
f    4   13
g    1   13
h    1   11
> duplicated(m)
[1] FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE  TRUE
> m[duplicated(m), ]
  [,1] [,2]
d    2   13
e    3   11
h    1   11

Bill Dunlap
TIBCO Software
wdunlap tibco.com



On Mon, Jul 21, 2014 at 7:54 AM, carol white  wrote:
> Hi,
> is it possible to find the duplicated rows of a matrix without a loop or i 
> have to loop over the rows? duplicated doesn't seem to be helpful
>
> Thanks
>
> Carol
>
>         [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] standard error of survfit.coxph()

[array chip]

Terry, I figured out that variance of log(-log(S)) should be (1/H^2)var(H), not 
(1/S^2)var(H)!

Thanks

John






e...@mayo.edu>; "r-help@r-project.org"  
Sent: Monday, July 21, 2014 11:41 AM
Subject: Re:  standard error of survfit.coxph()



Dear Terry,

I was trying to use your explanation of the standard error estimate from 
survfit.coxph() to verify the standard error estimates for the method of 
log(log(S)), but couldn't get the estimates correct. Here is an example using 
the lung dataset:

> fit<-coxph(Surv(time,status)~wt.loss,lung)
> surv<-survfit(fit,newdata=lung[2:5,])$surv

> logstd<-survfit(fit,newdata=lung[2:5,])$std.err  
> logstd.time<-survfit(fit,newdata=lung[2:5,])$time


    ## std error of cumulative hazard at time=60

> logstd<-logstd[logstd.time==60,]
> logstd
[1] 0.01965858 0.01965858 0.01941871 0.01969248

    ## survival S at months 60

> surv<-surv[logstd.time==60,]  
> surv
    2 3 4 5 
0.9249238 0.9249238 0.9253038 0.9263392

Since survfit()$std.err was for cumulative hazard H=log(S), thus based on your 
explanation below, the standard error for log-log method is: (1/S^2)var(H)

> loglogstd<-sqrt(1/surv^2*(logstd^2))
> loglogstd

 2  3  4  5 
0.02125427 0.02125427 0.02098631 0.02125839

    ## the upper confidence interval should be
> exp(-exp(log(-log(surv))-qnorm(0.975)*loglogstd))
    2
 3 4 5 
0.9278737 0.9278737 0.9282031 0.9292363

But this is different from the output using summary.survfit:

> summary(survfit(fit,newdata=lung[2:5,],conf.type='log-log'),times=60)$upper[1,]
    2 3 4 5 
0.9534814 0.9534814 0.9535646 0.9548478

Can you please suggest what I did wrong in the calculation?

Thanks very much!

John





To: "Therneau, Terry M., Ph.D." ; "r-help@r-project.org" 
 
Sent: Monday, June 30, 2014 10:46 AM
Subject: Re:  standard error of survfit.coxph()



[[elided Yahoo spam]]

John




 From: "Therneau, Terry M., Ph.D." 

Sent: Monday, June 30, 2014 6:04 AM
Subject: Re:  standard error of survfit.coxph()


1. The computations "behind the scenes" produce the variance of the cumulative 
hazard. 
This is true for both an ordinary Kaplan-Meier and a Cox model.  
Transformations to other 
scales are done using simple Taylor series.

   H = cumulative hazard = log(S);  S=survival
   var(H) = var(log(S))  = the starting point
   S = exp(log(S)), so  var(S) is approx [deriv of exp(x)]^2 * var(log(S)) = 
S^2 var(H)
   var(log(log(S)) is approx (1/S^2) var(H)

2. At the time it was written, summary.survfit was used only for
 printing out the survival 
curve at selected times, and the audience for the printout wanted std(S).   
True, that was 
20 years ago, but I don't recall anyone ever asking for summary to do anything 
else.  Your 
request is not a bad idea.
   Note however that the primary impact of using log(S) or S or log(log(S)) 
scale is is on 
the confidence intervals, and they do appear per request in the summary output.

Terry T.


On 06/28/2014 05:00 AM, r-help-requ...@r-project.org wrote:
> Message: 9
> Date: Fri, 27 Jun 2014 12:39:29 -0700

> To:"r-help@r-project.org"  
> Subject: [R] standard error of survfit.coxph()
> Message-ID:
>     <1403897969.91269.yahoomail...@web122906.mail.ne1.yahoo.com>
> Content-Type: text/plain
>
> Hi, can anyone help me to understand the standard errors printed in the 
> output of survfit.coxph()?
>
> time<-sample(1:15,100,replace=T)
>
> status<-as.numeric(runif(100,0,1)<0.2)
> x<-rnorm(100,10,2)
>
>
 fit<-coxph(Surv(time,status)~x)
> ??? ### method 1
>
> survfit(fit, newdata=data.frame(time=time,status=status,x=x)[1:5,], 
> conf.type='log')$std.err
>
>  [,1]??? [,2]??? [,3]???
 [,4]?? [,5]
> ?[1,] 0.0 0.0 0.0 0.0 0.
> ?[2,] 0.008627644 0.008567253 0.008773699 0.009354788 0.01481819
> ?[3,] 0.008627644 0.008567253 0.008773699 0.009354788 0.01481819
> ?[4,] 0.013800603 0.013767977 0.013889971 0.014379928 0.02353371
> ?[5,] 0.013800603 0.013767977 0.013889971 0.014379928 0.02353371
> ?[6,] 0.013800603 0.013767977 0.013889971 0.014379928 0.02353371
> ?[7,] 0.030226811 0.030423883 0.029806263 0.028918817 0.05191161
> ?[8,] 0.030226811 0.030423883 0.029806263 0.028918817 0.05191161
> ?[9,] 0.036852571 0.037159980 0.036186931 0.034645002 0.06485394
> [10,] 0.044181716 0.044621159 0.043221145 0.040872939 0.07931028
> [11,] 0.044181716 0.044621159 0.043221145 0.040872939 0.07931028
> [12,] 0.055452631 0.056018832 0.054236881 0.051586391 0.10800413
> [13,] 0.070665160 0.071363749 0.069208056 0.066655730
 0.14976433
> [14,] 0.124140400 0.125564637 0.121281571 0.118002021 0.30971860
> [15,] 0.173132357 0.175309455 0.168821266 0.164860523 0.46393111
>
> survfit(fit, newdata=data.frame(time=time,status=status,x=x)[1:5,], 
> conf.type='log')$time
> ?[1]? 1? 2? 3? 4? 5? 6?

Re: [R] Application design.

[Richard M. Heiberger]

Since your boss is Excel based, you might want to give the appearance
of staying in
that environment.  Take a look at RExcel and RWord, both at rcom.univie.ac.at
RExcel is a seamless integration of Excel and R.

See the book R through Excel that Erich Neuwirth (the author of
RExcel) and I wrote.
http://www.springer.com/mathematics/computational+science+%26+engineering/book/978-1-4419-0051-7

Rich

On Mon, Jul 21, 2014 at 10:24 PM, John McKown
 wrote:
> I'm designing an R based application for my boss. It's not much, but
> it might save him some time. What it will be doing is reading data
> from an MS-SQL database and creating a number of graphs. At present,
> he must log into one server to run a vendor application to display the
> data in a grid. He then cuts this data and pastes it into an Excel
> spreadsheet. He then generates some graphs in Excel. Which he then
> cuts and pastes into a Power Point presentation. Which is the end
> result for distribution to others up the food chain.
>
> What I would like to do is read the MS-SQL data base using RODBC and
> create the graphs using ggplot2 instead of using Excel. I may end up
> being told to create an Excel file as well.
>
> My real question is organizing the R programs to do this. Basically
> what I was thinking of was a "master" program. It does the ODBC work
> and fetches the data into one, or more, data.frames. I was then
> thinking that it would be better to have separate source files for
> each graph produced. I would use the source() function in the "master"
> R program to load & execute each one in order. Is this a decent
> origination? Or would it be better for each "create a graph" R file to
> really just define a unique function which the "master" program would
> then invoke? I guess this latter would be a good way to keep the
> workspace "clean" since all the variables in the functinon would "go
> away" when the function ended.
>
> I guess what I'm asking is how others organize the R applications. Oh,
> I plan for this to be run by my boss by double clicking on the
> "master" R source file, which I will associate with the Rscript
> program in Windows. Yes, this is Windows based .
>
> Appreciate your thoughts. Especially if I'm really off track.
>
> --
> There is nothing more pleasant than traveling and meeting new people!
> Genghis Khan
>
> Maranatha! <><
> John McKown
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Application design.

[John McKown]

I'm designing an R based application for my boss. It's not much, but
it might save him some time. What it will be doing is reading data
from an MS-SQL database and creating a number of graphs. At present,
he must log into one server to run a vendor application to display the
data in a grid. He then cuts this data and pastes it into an Excel
spreadsheet. He then generates some graphs in Excel. Which he then
cuts and pastes into a Power Point presentation. Which is the end
result for distribution to others up the food chain.

What I would like to do is read the MS-SQL data base using RODBC and
create the graphs using ggplot2 instead of using Excel. I may end up
being told to create an Excel file as well.

My real question is organizing the R programs to do this. Basically
what I was thinking of was a "master" program. It does the ODBC work
and fetches the data into one, or more, data.frames. I was then
thinking that it would be better to have separate source files for
each graph produced. I would use the source() function in the "master"
R program to load & execute each one in order. Is this a decent
origination? Or would it be better for each "create a graph" R file to
really just define a unique function which the "master" program would
then invoke? I guess this latter would be a good way to keep the
workspace "clean" since all the variables in the functinon would "go
away" when the function ended.

I guess what I'm asking is how others organize the R applications. Oh,
I plan for this to be run by my boss by double clicking on the
"master" R source file, which I will associate with the Rscript
program in Windows. Yes, this is Windows based .

Appreciate your thoughts. Especially if I'm really off track.

-- 
There is nothing more pleasant than traveling and meeting new people!
Genghis Khan

Maranatha! <><
John McKown

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Help with SEM package - model significance

[John Fox]

Dear Bernado,

 This isn't really a suitable topic to pursue on the r-help list, so I'll just 
comment briefly:

On Mon, 21 Jul 2014 17:34:52 -0700
 Bernardo Santos  wrote:
> Hi John,
> 
> Thanks for your reply (1 month later lol).
> In fact maybe the point is that I do not understand exactly the role of 
> latent variables (what they are, and how to define them in R) in SEM.
> Do you have any suggestion of easy basic literature on SEM that can help me 
> with that?
> Most things I have read are old (and use some different statistics programs) 
> and offer examples too far from ecology, than I had some difficulties to 
> understand the method in general.
> 

I have some materials at 
 from a 
recent workshop on SEMs, including some reading and suggestions for reading, 
but not I'm afraid from ecology.
 
> But, as far as I understood, SEM is like a simple multiple regression (linear 
> model), but that takes into account the relation of different variables 
> simultaneously, isn't that?
> 

In SEMs the response variable from one regression equation can be an 
explanatory variable in another, and the models can incorporate latent 
variables, which aren't measured directly, but rather indirectly through their 
observable effects ("indicators") or even in some cases through their 
observable causes.

I hope this helps,
 John

> Thank you very much.
> Best regards,
> 
> Bernardo
> 
> 
> 
> Em Segunda-feira, 16 de Junho de 2014 8:40, John Fox  
> escreveu:
>  
> 
> 
> Dear Bernardo,
> 
> The df for the LR chisquare over-identification test come not from the number 
> of observations, but from the difference between the number of observable 
> variances and covariances, on the one hand, and free parameters to estimate, 
> on the other. In your case, these numbers are equal, and so df = 0. The LR 
> chisquare for a just-identified model is also necessarily 0: the model 
> perfectly reproduces the covariational structure of the observed variables. 
> 
> R (and most statistical software) by default writes very small and very large 
> numbers in scientific format. In your case, -2.873188e-13 = -2.87*10^-13, 
> that is, 0 within rounding error. You can change the way numbers are printed 
> with the R scipen option.
> 
> Some other observations:
> 
> (1) Your model is recursive and has no latent variables; you would get the 
> same estimates from OLS regression using lm().
> 
> (2) For quite some time now, the sem package has included specifyEquations() 
> as a more convenient way of specifying a model, in preference to 
> specifyModel(). See ?specifyEquations.
> 
> (3) You don't have to specify the error variances directly; 
> specifyEquations(), or specifyModel(), will supply them.
> 
> I hope this helps,
> John
> 
> 
> 
> John Fox, Professor
> McMaster University
> Hamilton, Ontario, Canada
> http://socserv.mcmaster.ca/jfox/
>     
>     
> On Sun, 15 Jun 2014 20:15:31 -0700 (PDT)
> Bernardo Santos  wrote:
> > Dear all, 
> > 
> > I used "sem" function from the package SEM to fit a model. However, I 
> > cannot say if the model is correspondent to the data or not (chisquare 
> > test).
> > I used the commands:
> > 
> > model1 <- specifyModel()
> > estadio -> compflora, a1, NA
> > estadio -> compfauna, a2, NA
> > estadio -> interacoesobs, a3, NA
> > compflora -> compfauna, b1, NA
> > compflora -> interacoesobs, b2, NA
> > compfauna -> interacoesobs, c1, NA
> > estadio <-> estadio, e1, NA
> > compflora <-> compflora, e2, NA
> > compfauna <-> compfauna, e3, NA
> > interacoesobs <-> interacoesobs, e4, NA
> > 
> > sem1 <- sem(model1, cov.matrix, length(samples))
> > summary(sem1)
> > 
> > and I got the result:
> > 
> > Model Chisquare =  -2.873188e-13   Df =  0 Pr(>Chisq) = NA AIC =  20 BIC =  
> > -2.873188e-13 Normalized Residuals Min.   1st Qu.    Median      Mean   3rd 
> > Qu.      Max. 
> > 0.000e+00 0.000e+00 2.957e-16 3.193e-16 5.044e-16 8.141e-16  R-square for 
> > Endogenous Variables compflora     compfauna interacoesobs  0.0657        
> > 0.1056        0.2319  Parameter Estimates Estimate     Std Error    z value 
> >    Pr(>|z|)                                    
> > a1 3.027344e-01 1.665395e-01 1.81779316 6.909575e-02 compflora <--- estadio 
> >          
> > a2 2.189427e-01 1.767404e-01 1.23878105 2.154266e-01 compfauna <--- estadio 
> >          
> > a3 7.314192e-03 1.063613e-01 0.06876742 9.451748e-01 interacoesobs <--- 
> > estadio      
> > b1 2.422906e-01 1.496290e-01 1.61927587 1.053879e-01 compfauna <--- 
> > compflora        
> > b2 3.029933e-01 9.104901e-02 3.32780446 8.753328e-04 interacoesobs <--- 
> > compflora    
> > c1 4.863368e-02 8.638177e-02 0.56300857 5.734290e-01 interacoesobs <--- 
> > compfauna    
> > e1 6.918133e+04 1.427102e+04 4.84767986 1.249138e-06 estadio <--> estadio   
> >          
> > e2 9.018230e+04 1.860319e+04 4.84767986 1.249138e-06 compflora <--> 
> > com

Re: [R] Semi Markov warnings ( for dummies)

[Rolf Turner]

If you are going to drive a car you should learn to drive.  Read the 
basic intro material for R.


Where does the function "readdata" come from?  The syntax you use 
(having an assignment inside the function call) is let us say unorthodox 
and highly inadvisable.


Since you got an error reading in your data, where did the data that you 
are trying to analyse come from?


Don't call your data "data".  See fortune(77).

You don't say where the function semiMarkov() comes from; presumably 
from the package "semiMarkov".  Say so.  Do not expect your readers to 
be telepathic.


Don't use Excel to do *anything* (except store data, if you *must*).

I certainly can't help you to debug your errors/warnings without seeing 
your data.  I doubt that anyone else can either.  (I *could* be wrong 
about that; there are some *very* clever and intuitive people on this list.)


cheers,

Rolf Turner

On 22/07/14 02:04, M.A. Pet wrote:

*Hello,*

I never worked with R before my supervisor asks me to run a semiMarkov
analysis a month ago. After a long struggle, to date, the code works, but I
still get some warnings. However, because of my lack of knowledge in R I am
not possible to figure out the problems or say anything about the influence
of these warnings on my outcome. Hopefully, someone would help me with
these ( I think basic) questions.


What is the case? I want to do a semiMarkov analysis with 3 states (state
n, state s and state e). Wherefore I want to run the analysis to see
whether there is a difference in Hazard Ratio between transitions nà e and sà e.
Therefore, I’ve got data of nearly 60 persons. In excel, every worksheet
reflects a person.  To answer my research question I tried to run the
script for one person.

Unfortunately, I get the following error:


# read data and convert to correct format (and print tail).



data <- readdata( input <- 'E1_tijd_stress.csv' )


Error in `$<-.data.frame`(`*tmp*`, "state", value = "s") :

   replacement has 1 row, data has 0

Called from: `$<-`(`*tmp*`, "state", value = "s")



Furthermore, I used a macro in excel to merge these worksheets into one
worksheet ( all persons in one file) . However, after running the script I
saw the following warnings:


# semi-Markov model without covariates.



fit.semi.markov <- semiMarkov( data = data, states = states.semi.markov, mtrans 
= mtrans.semi.markov )




Iter: 1 fn: 11224.5658  Pars:  1000.0 412.24725 261.14293
66.30789   1.67507   5.59420   0.55746   1.37715   0.59734   0.48604
0.73281   0.74205

Iter: 2 fn: 11224.5658  Pars:  1000.0 412.24425 261.14508
66.30635   1.67506   5.59420   0.55746   1.37716   0.59734   0.48604
0.73281   0.74206

solnp--> Completed in 2 iterations

Warning messages:

1: In .safefunx(tmpv, .solnp_fun, .env, ...) :

solnp-->warning: Inf detected in function call...check your function


These repeated for 7 times. Or the following warning:


fit.semi.markov <- semiMarkov( data = data, states = states.semi.markov, mtrans 
= mtrans.semi.markov )


Iter: 1 fn: 444562.8860 Pars:  1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0 0.47882 0.72691 0.74627

solnp--> Completed in 1 iterations

Error in svd(X) : infinite or missing values in 'x'

In addition: There were 29 warnings (use warnings() to see them)

29: In .safefunx(tmpv, .solnp_fun, .env, ...) :

solnp-->warning: NaN detected in function call...check your function


Furthermore, I want to control the analysis for some covariates.


# semi-Markov model with covariates



fit2 <- semiMarkov( data = data,cov= as.data.frame(data$stress.score), states = 
states, mtrans = mtrans )


Error in mtrans[as.numeric(substring(trans.hh[i], first = 1, last = 1)),  :

   subscript out of bounds


Is someone able to tell me what these warnings exactly mean? Or someone I
can contact to discuss about this? I understand inf and NaN for example but
I am not able to figure out what to change in my data format of script (
the problem is, some of my data run without problems, some does not).
Hopefully I give enough information, otherwise let me know.

Thank you!

[[alternative HTML version deleted]]



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.




--
Rolf Turner
Technical Editor ANZJS

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] roxygen2

[Brian Diggs]

On 7/20/2014 12:50 AM, Kevin Kunzmann wrote:

Hi,

I have developed a package and would like to switch documentation to
roxygen2 from manual :) However

 >roxygen2::roxygenize()
 First time using roxygen2 4.0. Upgrading automatically...
 Loading required package: nleqnslv
 Error in file(filename, "r", encoding = encoding) :
cannot open the connection
 In addition:Warning messages:
 1: In library(package, lib.loc = lib.loc, character.only = TRUE, 
logical.return = TRUE,  :
there is no package called 'nleqnslv'
 2: In file(filename, "r", encoding = encoding) :
cannot open file 'R/misc.R': No such file or directory

gives me this, any ideas? Running R 3.1.1 and the file R/misc.R DOES
exist. working directory is package directory.


I'd focus on the other warning message first: "there is no package 
called 'nleqnslv'". I don't see a package by that name on CRAN, though 
google helpfully suggested 'nleqslv' as something that does exist.


Do you have a line `library("nleqnslv")` in `R/misc.R`? If so, that may 
be the problem (also, if you do, you should not in general have 
`library` calls in package code; package interdependencies are handled 
by the DEPENDS and IMPORTS fields of the DESCRIPTION file). Or do you 
have a reference to "nleqnslv" in the DESCRIPTION file?


I would get the first warning message cleared up first (especially since 
the second one doesn't make sense to you) and see if the other one goes 
away.



Best,

Kevin


--
Brian S. Diggs, PhD
Senior Research Associate, Department of Surgery
Oregon Health & Science University

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] Include plotting symbals pch 16 and 17 into captions / text in pdf graph

[David Winsemius]

It's a hack but this works:

cairo_pdf("utftext.pdf", family="Calibri")
plot(0:1,0:1, type="n")
points(0.2, 1, pch=16)
mtext(
   text = "pch=16 (\U25CF)",
   side = 3,
   at   = 0.2,
   line = 1,
   )
points(0.8, 1, pch=17); points(0.85, 1.1, pch=17, xpd=TRUE)
# par(xpd=TRUE) lets plotting go outside the usual clipping box.

mtext(
   text = "pch=17 (  )",  # makes room for the triangle
   side = 3,
   at   = 0.8,
   line = 1
   )
dev.off()

-- 
David.

On Jul 21, 2014, at 5:10 AM, Rainer M Krug wrote:

> 
> Hi
> 
> I want to include the plotting symbols pch=16 and pch=17 in text in
> the graph used as labels.
> 
> This works:
> 
> --8<---cut here---start->8---
> plot(0:1,0:1, type="n")
> points(0.2, 1, pch=16)
> mtext(
>text = "pch=16 (\U25CF)",
>side = 3,
>at   = 0.2,
>line = 1,
>)
> points(0.8, 1, pch=17)
> mtext(
>text = "pch=17 (\U25B2)",
>side = 3,
>at   = 0.8,
>line = 1
>)
> --8<---cut here---end--->8---
> 
> But I would need this as a pdf and in Calibri font, which does not
> include the UTF symbol \U25B2 (the triangle). I got that far:
> 
> --8<---cut here---start->8---
> 
> cairo_pdf("utftext.pdf", family="Calibri")
> plot(0:1,0:1, type="n")
> points(0.2, 1, pch=16)
> mtext(
>text = "pch=16 (\U25CF)",
>side = 3,
>at   = 0.2,
>line = 1,
>)
> points(0.8, 1, pch=17)
> mtext(
>text = "pch=17 (\U25B2)",
>side = 3,
>at   = 0.8,
>line = 1
>)
> dev.off()
> --8<---cut here---end--->8---
> 
> which results in the attached pdf (hope it comes through). If not: the
> utf symbol for the upward error is displayed as an empty square with a
> question mark in it).
> 
> Is there a way that I can show the plotting symbol pch=17 in the caption
> when using this font?
> 
> Thanks,
> 
> Rainer
> 
> -- 
> Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology, 
> UCT), Dipl. Phys. (Germany)
> 
> Centre of Excellence for Invasion Biology
> Stellenbosch University
> South Africa
> 
> Tel :   +33 - (0)9 53 10 27 44
> Cell:   +33 - (0)6 85 62 59 98
> Fax :   +33 - (0)9 58 10 27 44
> 
> Fax (D):+49 - (0)3 21 21 25 22 44
> 
> email:  rai...@krugs.de
> 
> Skype:  RMkrug
> 
> PGP: 0x0F52F982
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

David Winsemius
Alameda, CA, USA

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] duplicated rows of a matrix

[William Dunlap]

duplicated(x), for vector or matrix x, flags any value (row for
matrices) previously seen in x.  To flag all duplicated values (rows
for matrices) you can use the following allDups() function.

   allDups <- function(x) duplicated(x) | duplicated(x, fromLast=TRUE)

> In your example, rows 2,3,4,5,8

Don't you want row 1 also, since rows 1 and 8 contain c(1,11)?  Or I
am again misinterpreting what you want?

Bill Dunlap
TIBCO Software
wdunlap tibco.com


On Mon, Jul 21, 2014 at 12:51 PM, carol white  wrote:
> I need that duplicated indicate all row indices (occurences) that are
> duplicated. In your example, rows 2,3,4,5,8
>
> Thanks.
>
> Carol
>
>
> On Monday, July 21, 2014 9:17 PM, William Dunlap  wrote:
>
>
> Can you give an example of duplicated() not working on the rows of a matrix?
>
> Here is an example where it does work:
>> m <- cbind(c(a=1,b=2,c=3,d=2,e=3,f=4,g=1,h=1), c(11,13,11,13,11,13,13,11))
>> class(m)
> [1] "matrix"
>> m
>   [,1] [,2]
> a1  11
> b2  13
> c3  11
> d2  13
> e3  11
> f4  13
> g1  13
> h1  11
>> duplicated(m)
> [1] FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE  TRUE
>> m[duplicated(m), ]
>   [,1] [,2]
> d2  13
> e3  11
> h1  11
>
> Bill Dunlap
> TIBCO Software
> wdunlap tibco.com
>
>
> On Mon, Jul 21, 2014 at 7:54 AM, carol white  wrote:
>> Hi,
>> is it possible to find the duplicated rows of a matrix without a loop or i
>> have to loop over the rows? duplicated doesn't seem to be helpful
>>
>> Thanks
>>
>> Carol
>
>>
>>[[alternative HTML version deleted]]
>>
>> __
>> R-help@r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide
>> http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>
>
>

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] anova.lme

[Robert Lynch]

I would like to know the sum of squares for each term in my model. I used
the following call to fit the model

fit.courseCross <- lme(fixed= zGrade ~ Rep  + ISE
+P7APrior+Female+White+HSGPA+MATH+Years+Course+Course*P7APrior ,
  random= ~1|SID,
  data = Master.complete[Master.complete$Course != "P7A",])

and called an anova on it and get:

anova(fit.courseCross)
numDF denDF   F-value p-value
(Intercept) 1 58161 1559.6968  <.0001
Rep 1 58161  520.7263  <.0001
ISE 1  6266   21.3713  <.0001
P7APrior2 58161  358.4827  <.0001
Female  1  6266   89.2614  <.0001
White   1  6266  235.9984  <.0001
HSGPA   1  6266 1156.4116  <.0001
MATH1  6266 1036.1354  <.0001
Years   1 58161  407.6096  <.0001
Course 12 58161   68.9875  <.0001
P7APrior:Course24 58161   10.2464  <.0001

The documentation for anova.lme says:

When only one fitted model object is present, a data frame with the sums of
squares, numerator degrees of freedom, denominator degrees of freedom,
F-values, and P-values for Wald tests for the terms in the model (when
Terms and L are NULL), a combination of model terms (when Terms in not
NULL), or linear combinations of the model coefficients (when L is not
NULL).

noticeably absent is the sum of squares.

How do I get them?
Robert

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] duplicated rows of a matrix

[William Dunlap]

Can you give an example of duplicated() not working on the rows of a matrix?

Here is an example where it does work:
> m <- cbind(c(a=1,b=2,c=3,d=2,e=3,f=4,g=1,h=1), c(11,13,11,13,11,13,13,11))
> class(m)
[1] "matrix"
> m
  [,1] [,2]
a1   11
b2   13
c3   11
d2   13
e3   11
f4   13
g1   13
h1   11
> duplicated(m)
[1] FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE  TRUE
> m[duplicated(m), ]
  [,1] [,2]
d2   13
e3   11
h1   11

Bill Dunlap
TIBCO Software
wdunlap tibco.com


On Mon, Jul 21, 2014 at 7:54 AM, carol white  wrote:
> Hi,
> is it possible to find the duplicated rows of a matrix without a loop or i 
> have to loop over the rows? duplicated doesn't seem to be helpful
>
> Thanks
>
> Carol
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] standard error of survfit.coxph()

[array chip]

Dear Terry,

I was trying to use your explanation of the standard error estimate from 
survfit.coxph() to verify the standard error estimates for the method of 
log(log(S)), but couldn't get the estimates correct. Here is an example using 
the lung dataset:

> fit<-coxph(Surv(time,status)~wt.loss,lung)
> surv<-survfit(fit,newdata=lung[2:5,])$surv

> logstd<-survfit(fit,newdata=lung[2:5,])$std.err  
> logstd.time<-survfit(fit,newdata=lung[2:5,])$time


    ## std error of cumulative hazard at time=60

> logstd<-logstd[logstd.time==60,]
> logstd
[1] 0.01965858 0.01965858 0.01941871 0.01969248

    ## survival S at months 60

> surv<-surv[logstd.time==60,]  
> surv
    2 3 4 5 
0.9249238 0.9249238 0.9253038 0.9263392

Since survfit()$std.err was for cumulative hazard H=log(S), thus based on your 
explanation below, the standard error for log-log method is: (1/S^2)var(H)

> loglogstd<-sqrt(1/surv^2*(logstd^2))
> loglogstd

 2  3  4  5 
0.02125427 0.02125427 0.02098631 0.02125839

    ## the upper confidence interval should be
> exp(-exp(log(-log(surv))-qnorm(0.975)*loglogstd))
    2 3 4 5 
0.9278737 0.9278737 0.9282031 0.9292363

But this is different from the output using summary.survfit:

> summary(survfit(fit,newdata=lung[2:5,],conf.type='log-log'),times=60)$upper[1,]
    2 3 4 5 
0.9534814 0.9534814 0.9535646 0.9548478

Can you please suggest what I did wrong in the calculation?

Thanks very much!

John





To: "Therneau, Terry M., Ph.D." ; "r-help@r-project.org" 
 
Sent: Monday, June 30, 2014 10:46 AM
Subject: Re:  standard error of survfit.coxph()



[[elided Yahoo spam]]

John




 From: "Therneau, Terry M., Ph.D." 

Sent: Monday, June 30, 2014 6:04 AM
Subject: Re:  standard error of survfit.coxph()


1. The computations "behind the scenes" produce the variance of the cumulative 
hazard. 
This is true for both an ordinary Kaplan-Meier and a Cox model.  
Transformations to other 
scales are done using simple Taylor series.

   H = cumulative hazard = log(S);  S=survival
   var(H) = var(log(S))  = the starting point
   S = exp(log(S)), so  var(S) is approx [deriv of exp(x)]^2 * var(log(S)) = 
S^2 var(H)
   var(log(log(S)) is approx (1/S^2) var(H)

2. At the time it was written, summary.survfit was used only for
 printing out the survival 
curve at selected times, and the audience for the printout wanted std(S).   
True, that was 
20 years ago, but I don't recall anyone ever asking for summary to do anything 
else.  Your 
request is not a bad idea.
   Note however that the primary impact of using log(S) or S or log(log(S)) 
scale is is on 
the confidence intervals, and they do appear per request in the summary output.

Terry T.


On 06/28/2014 05:00 AM, r-help-requ...@r-project.org wrote:
> Message: 9
> Date: Fri, 27 Jun 2014 12:39:29 -0700

> To:"r-help@r-project.org"  
> Subject: [R] standard error of survfit.coxph()
> Message-ID:
>     <1403897969.91269.yahoomail...@web122906.mail.ne1.yahoo.com>
> Content-Type: text/plain
>
> Hi, can anyone help me to understand the standard errors printed in the 
> output of survfit.coxph()?
>
> time<-sample(1:15,100,replace=T)
>
> status<-as.numeric(runif(100,0,1)<0.2)
> x<-rnorm(100,10,2)
>
> fit<-coxph(Surv(time,status)~x)
> ??? ### method 1
>
> survfit(fit, newdata=data.frame(time=time,status=status,x=x)[1:5,], 
> conf.type='log')$std.err
>
>  [,1]??? [,2]??? [,3]???
 [,4]?? [,5]
> ?[1,] 0.0 0.0 0.0 0.0 0.
> ?[2,] 0.008627644 0.008567253 0.008773699 0.009354788 0.01481819
> ?[3,] 0.008627644 0.008567253 0.008773699 0.009354788 0.01481819
> ?[4,] 0.013800603 0.013767977 0.013889971 0.014379928 0.02353371
> ?[5,] 0.013800603 0.013767977 0.013889971 0.014379928 0.02353371
> ?[6,] 0.013800603 0.013767977 0.013889971 0.014379928 0.02353371
> ?[7,] 0.030226811 0.030423883 0.029806263 0.028918817 0.05191161
> ?[8,] 0.030226811 0.030423883 0.029806263 0.028918817 0.05191161
> ?[9,] 0.036852571 0.037159980 0.036186931 0.034645002 0.06485394
> [10,] 0.044181716 0.044621159 0.043221145 0.040872939 0.07931028
> [11,] 0.044181716 0.044621159 0.043221145 0.040872939 0.07931028
> [12,] 0.055452631 0.056018832 0.054236881 0.051586391 0.10800413
> [13,] 0.070665160 0.071363749 0.069208056 0.066655730
 0.14976433
> [14,] 0.124140400 0.125564637 0.121281571 0.118002021 0.30971860
> [15,] 0.173132357 0.175309455 0.168821266 0.164860523 0.46393111
>
> survfit(fit, newdata=data.frame(time=time,status=status,x=x)[1:5,], 
> conf.type='log')$time
> ?[1]? 1? 2? 3? 4? 5? 6? 7? 8? 9 10 11 12 13 14 15
>
> ??? ### method 2
>
> summary(survfit(fit, newdata=data.frame(time=time,status=status,x=x)[1:5,], 
> conf.type='log'),time=10)$std.err
>
> ? 1? 2? 3? 4? 5
> [1,] 0.04061384 0.04106186 0.03963184 0.03715246

[R] Error message for corAR1()

[Wilson, Jenny]

Hi,


I am trying to answer the see if density.km (response) is affected by Direction 
(continuous, integer), Layer (nominal with 12 levels) and direction (nominal 
with 8 levels). There is an interaction between Layer and Direction. 
Platform.field is a list of 9 different platforms and is being treated as a 
random effect.


I had previously ran this model without the correlation argument and checked 
the residuals in an acf plot which showed a high level of autocorrelation.


I am having issues applying the correlation argument into my model but keep 
getting an error message and am not sure what else to try. I should mention 
that Direction is being applied as the time covariate because each distance has 
an associated time stamp.


G2<-gamm(density.km~f.Layer+direction+s(Distance.A,by=f.Layer),
random=list(Platform.field=~1),corr=corAR1(form=~Distance.A|Platform.field/f.Layer),
family=poisson,data=fish1)


 Maximum number of PQL iterations:  20
iteration 1
Error in Initialize.corAR1(X[[2L]], ...) :
  covariate must have unique values within groups for "corAR1" objects?


Any help would be appreciated.


Thanks



[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] deviance as a goodness of fit in GLM

[Samantha PameLa]

Good day everybody, 
 
I'm a marine biologist student, working on my bachelor thesis and I'm stucked 
with a statistical doubt in the process, I hope someone here could help me. My 
thesis aims to understand which biological and environmental factors influences 
the male aggressive rate of male California sea lions. For that purpose I'm 
using GLM’s where the response variable is the male aggressive rate. Right now 
I am testing the goodness of fit of the global models and for that I'm using 
the deviances as a goodness of fit test. I calculated pseudoR2 (Zuur, 2009) in 
order to know the percentage of explanation of each candidate model. However 
I’m not sure how to choose the “good models”; since I am not sure over which 
percent of explanation indicates a “model with good fit”. For my data I am 
working with three different scenarios, and it seems that 20%, is a good value 
to could indicate the best models, but I’m not sure how to choose the value.
 
I thank you in advance for your time and the help you can give me.
 
Best regards,
 
Samantha.




  
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] duplicated rows of a matrix

[carol white]

Hi,
is it possible to find the duplicated rows of a matrix without a loop or i have 
to loop over the rows? duplicated doesn't seem to be helpful

Thanks

Carol

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Estimation of Zero Inflated Over dispersed Beta Binomial Using glamADMB()

[Rajibul Mian]

Dear All,

I have been facing problem running the following code by using
---glamADMB()--

glmmadmb(y_zibb~x+factor(z)+g, data= data_mis_model, family =
"betabinomial", link = "logit", zeroInflation=T)

where "y_zibb" contains zero inflated Beta Binomial response ,
"x" is a normal random variate
"z" is a binomial random variate
"g" is a exponential random variate

the error message

Error in glmmadmb(y_zibb ~ x + factor(z) + g, data = data_mis, family =
"betabinomial",  :
  The function maximizer failed (couldn't find STD file) Troubleshooting
steps include (1) run with 'save.dir' set and inspect output files; (2)
change run parameters: see '?admbControl'
In addition: Warning message:
running command 'C:\Windows\system32\cmd.exe /c
"C:/Users/rajibulmian/Documents/R/win-library/3.1/glmmADMB/bin/windows64/glmmadmb.exe"
-maxfn 500 -maxph 5 -noinit -shess' had status 22


I have tried the  "change run parameters: see '?admbControl'" in different
combinations but couldn't help.  I am giving part of the data with this
mail.

Any related suggestions would help. Thank you.

Best.
Rajibul Mian
Grad Student, Maths & STATS,
UWindsor, Canada.

> data_1.1
y_zibb   x zg
10 -0.70831528 1  1.602785996
23 -0.32331815 1  0.831407214
30 -2.05941215 1  0.650309590
40  0.22345640 1  2.013147059
54 -0.02260530 1  1.010770937
60  0.41063599 1  9.841200077
70 -0.27111008 1  4.775681308
80 -0.90399357 1  0.678469073
90 -0.21160307 1  0.587509666
10   0  0.51008466 1 13.319207785
-
-
-

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Generating nonlinear Poisson time series data

[Kristynn Sullivan]

We are attempting to create a short Poisson time series (between 10 and 50
datapoints) for a simulation. We want these time series to have no counts
of over 100 and not be zero-inflated. We also are trying generate various
nonlinearities, particularly of a cyclic nature. We have been attempting to
generate this data using higher order polynomials. An example with a
seventh-order trend (and a treatment effect):

Time <- 0:(T-1) ##T is the desired number of time points in the time series.
beta <- c(B0 = 1.585, B1 = Btrt, B2 = 1.229, B3 = -3.364e-01, B4 =
-6.610e-02, B5 = 2.697e-02, B6 = -2.905e-03, B7 =  1.304e-04, B8 =
-2.130e-06)
pmat <- cbind(const = 1, tx = tx, Time = Time, Time2 = Time^2, Time3 =
Time^3, Time4 = Time^4, Time5 = Time^5,Time6 = Time^6, Time7 = Time^7)
##Btrt is the treatment effect
 y <- pmat %*% beta
 y <- rpois(T, exp(y))

This code works. However, when manipulating the factors such as the length
of the time series, the same beta coefficients do not always produce the
same desirable properties - the most common problem being the end of the
time series results in a line of zeros. We are hoping to use a consistent
data generating method across simulation cells (e.g. not changing the set
of beta coefficients for each cell), so they are comparable in the analysis
stage of the simulation.

Does anyone have any other ways they create cyclic Poisson time series
data? Or good references for doing so?

Thanks in advance,

Kristynn Sullivan

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Weight, weight - do tell me

[Lenth, Russell V]

This is a question only about terminology.

Suppose I have data categorized by three factors A, B, and C, with cell means 
ybar_ijk and cell frequencies n_ijk, where I, j, and k index A, B, and C 
respectively. And suppose I want to summarize the results for factor A by 
computing some sort of weighted means WM_i, averaging over indices j and k with 
weights w_jk. Consider these four weighting schemes:

1. Use equal weights, w_jk = 1
2. Use weights of w_jk = n_+jk   (where "+" shows I summed over that index)
3. Use weights w_jk = n_+j+ * n_++k   (outer product of the one-factor 
marginal frequencies)
4. Use weights w_ijk = n_ijk   (only one where we use a different set of 
weights for each i)

Scheme 1 yields the "unweighted" or "least-squares" means, and scheme 4 yields 
the ordinary means for A, ignoring B and C altogether. Scheme 3 yields weighted 
averages over k of weighted averages over j (or vice versa).

My question is what to call these schemes, e.g., as a character argument in an 
R function. Preliminarily, I am calling them "equal", "proportional", "outer", 
and "actual". The first one is pretty obvious. But maybe there is some existing 
standard terminology for some or all of the others that I am unaware of (or 
have forgotten). Any suggestions?

Russ

Russell V. Lenth  -  Professor Emeritus
Department of Statistics and Actuarial Science   
The University of Iowa  -  Iowa City, IA 52242  USA   
Voice (319)335-0712 (Dept. office)  -  FAX (319)335-3017

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] help with column substaction with a twist

[Greg Snow]

Here is another approach in R (blatantly stealing Jim Holtman's code
to generate sample data):

> set.seed(1)
> n <- 100
> test <- data.frame(p = sample(10, n, TRUE)
+ , b = sample(10, n, TRUE)
+ )
> test$e <- sample(5, n, TRUE) + test$b  # make sure e > b
>
> tmp1 <- test$b - test$p
> tmp2 <- test$p - test$e
>
> test$dist <- pmax( tmp1, tmp2, 0 ) * sign( -tmp1 )
>
> head(test, 10)
p  b  e dist
1   3  7  9   -4
2   4  4  60
3   6  3  60
4  10 10 120
5   3  7  8   -4
6   9  3  63
7  10  2  55
8   7  5  61
9   7 10 12   -3
10  1  6 10   -5
>

You could also skip the 2 temporary variables and just code the
differences inside the pmax and sign functions (using tmp1 will mean
not having to do the same subtraction twice).  If you are happy with
the absolute difference then you can drop the "* sign( -tmp1 )" part.

This works because if p is less than b then tmp1 will be positive and
tmp2 will be negative.  If p is between b and e then both will be
negative (and therefore 0 will be greater than both). If p is greater
than e then tmp2 will be positive and tmp1 negative.  So the maximum
value (pmax computes this for each row/pair/triplet) will be the one
of interest.


On Sun, Jul 20, 2014 at 5:33 AM, Ubernerdy  wrote:
> Hello guys!
>
> I have been messing around with R for a while now, but this situation has
> me a bit stumped. I was unable to solve it by reading documentation.
>
> So I have this table (currently in Excel - could export it as csv) with
> values in 3 columns. Let's call them value P (for position), value B
> (beginning) and E (end). P represents the position of a mutation in the
> genome, B and E are the beginnings and ends of a nearby gene that either
> contains the mutation or not.
>
> I am trying to compute the distance between the mutation and the gene.
>
> If the mutation is contained in the gene, that is value P is greater than B
> and lesser than E, the result is 0.
>
> If the mutation is "left" of the gene, the distance is negative and is
> equal to P-B.
>
> If the mutation is "right" of the gene, the distance is positive and is
> equal to P-E.
>
> How would i achieve this in R?
>
> Regards and thanks, S.
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.



-- 
Gregory (Greg) L. Snow Ph.D.
538...@gmail.com

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] FW: two questions - function help and 32vs64 bit sessions

[Fowler, Mark]

-Original Message-
From: Fowler, Mark 
Sent: July 21, 2014 1:56 PM
To: 'Duncan Murdoch'
Subject: RE: [R] two questions - function help and 32vs64 bit sessions

The server doesn't shut down, it just kind of shuts everybody else down.
But your mention of TEMPDIR jostled some old memory cells awake, and
essentially answered the question. I do [did] a similar thing with
temporary directories as you discuss, just another directory location.
But that location (R_USER) disappeared when I upgraded (one or all of R
to 2 to 3, OS XP to 7, machine 34 to 62 bit). So R has to create a temp
dir by default. I think that default directory gets derailed during
network-wide system scans, thus Tuesdays and weekends for me.

Sys.getenv('R_USER')
[1] "C:\\Users\\fowlerm\\Documents"
> tempdir()
[1] "C:\\Users\\fowlerm\\AppData\\Local\\Temp\\1\\RtmpeqbLeR"

-Original Message-
From: Duncan Murdoch [mailto:murdoch.dun...@gmail.com]
Sent: July 21, 2014 11:55 AM
To: Fowler, Mark; r-h...@stat.math.ethz.ch
Subject: Re: [R] two questions - function help and 32vs64 bit sessions

On 21/07/2014 9:40 AM, Fowler, Mark wrote:
> Hi Duncan,
> I tried your suggestion, but no luck. The first error is no surprise, 
> it just confirms the address is lost. The second line suggests it 
> worked, but it didn't. The session is still remembering the original
address.
>
> > tools::startDynamicHelp(FALSE) # shut it down
> Warning message:
> In file(out, "wt") :
>cannot open file
> 'C:\Users\fowlerm\AppData\Local\Temp\1\RtmpK09I4B\Rhttpd26c04742884': 
> No such file or directory
> > tools::startDynamicHelp(TRUE)  # start it up
> starting httpd help server ... done

This is a different error than I thought you described.  I thought
something had shut down the server, but it looks like something has 
deleted your temporary directory.   You might be able to tell your OS 
not to do that (if it was the OS that did), or you can put your
temporary directory in a location that is less likely to get deleted,
e.g. by setting the TMPDIR environment variable to somewhere else before
you start R.  For example, I typically run with TMPDIR set to C:/temp
when I'm debugging things.

Duncan Murdoch

>
> -Original Message-
> From: Duncan Murdoch [mailto:murdoch.dun...@gmail.com]
> Sent: July 16, 2014 10:47 AM
> To: Fowler, Mark; r-h...@stat.math.ethz.ch
> Subject: Re: [R] two questions - function help and 32vs64 bit sessions
>
> On 14/07/2014, 11:42 AM, Fowler, Mark wrote:
> > Hello,
> >
> >
> >
> > Two unrelated questions, and neither urgent.
> >
> >
> >
> > Windows 7, R 3.0.1. Using R Console, no fancy interface.
> >
> >
> >
> > The function help ultimately becomes lost to a session kept running 
> > for extended periods (days). I.e. with a new session if you invoke 
> > the
>
> > Help menu 'R functions (txt)...' it activates the html help and goes

> > to the named function page. This will work fine for at least a day, 
> > but typically the next day invoking the help menu in this fashion 
> > will
>
> > fail, as R looks for a temporary address it creates on your
computer.
> > This gets lost, possibly due to network administration activity. So 
> > then I save and start another session with same Rdata. Trivial 
> > enough but irritating. Anybody know how to restore the 'link'
> > without ending and restarting the session?
>
> I never have sessions that last that long, so I haven't tried this, 
> but I'd expect you could restart the help system in this way:
>
> tools::startDynamicHelp(FALSE) # shut it down
> tools::startDynamicHelp(TRUE)  # start it up
>
> Duncan Murdoch
>
> >
> >
> >
> > I have a mix of 32-bit and 64-bit requirements, with 64 the default.

> > I
>
> > became used to starting R sessions directly from the appropriate 
> > Rdata
>
> > workspaces. With the latest version I need to start from the generic

> > icon and then load the workspace if I want 32-bit. Anybody know a 
> > way to make the Rdata files keep track of which bit version they 
> > work with, or some trick that accomplishes the same objective? The 
> > 32-bit requirement is usually just RODBC, and the need for it is 
> > scattered over lots of workspaces. Again, trivial but a nuisance. A 
> > more pragmatic motive is to not oblige users of applications to 
> > think about
>
> > it. Any way to make a session switch 'bits' with a source file?
> >
> >
> >
> > Mark Fowler
> > Population Ecology Division
> > Bedford Inst of Oceanography
> > Dept Fisheries & Oceans
> > Dartmouth NS Canada
> > B2Y 4A2
> > Tel. (902) 426-3529
> > Fax (902) 426-9710
> > Email mark.fow...@dfo-mpo.gc.ca 
> >
> >
> >
> >
> >
> >
> > [[alternative HTML version deleted]]
> >
> > __
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> > http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducibl

Re: [R] Survival Analysis with an Historical Control

[Andrews, Chris]

Hi Paul,
Sorry for the delayed reply.  I was away last week. I'm not clear what you want 
confirmed about your approach.

(a) "20.57"- computing the rejection region of the analysis.  The formulas 
implemented at the addresses you gave in your original post are from a 
reputable source - Lawless (1982) - at least that is claimed in the help file 
(http://www.swogstat.org/stat/Public/Help/survival1.html).  It seems that 
another person derived 20.57 from some combination of input.  I don't see how 
to back calculate what the input was from the information provided.  Perhaps 
elsewhere in your study protocol there is discussion of accrual time, et al. 
that can help you.

(b) "16" - the value of the parameter in the null hypothesis is very important 
as you noted in your response to Terry.  But this is not really a statistical 
question. It may be derived from historical data, expert opinion, regulatory 
mandate, or some combination of these and other factors.  Presumably this study 
was undertaken because 16 was an important number to somebody.

(c) performing the one-sample survival analysis itself.  This is what I did 
with the data you provided.

# Non-parametric
(km <- survfit(Surv(WKS, 1-CENS) ~ 1, data=hsv, type="kaplan-meier", 
conf.type="log", conf.int=0.9))
summary(km)

# Compare to median survival = 16
# (Used 90% CI above to get 0.05 one sided test here)
quantile(km, prob=0.5)$lower > 16

# Parametric
(paraExp <- survreg(Surv(WKS, 1-CENS) ~ 1, data=hsv, dist="exponential"))
summary(paraExp)

# Compare to median survival = 16
# That is, compare to beta0 = log(16/log(2)) = 3.1391...
# one sided test
pnorm(c((coef(paraExp) - log(16/log(2))) / sqrt(vcov(paraExp))), 
lower.tail=FALSE)


Chris

-Original Message-
From: Paul Miller [mailto:pjmiller...@yahoo.com] 
Sent: Thursday, July 10, 2014 8:59 AM
To: Andrews, Chris
Cc: r-help@r-project.org
Subject: RE: [R] Survival Analysis with an Historical Control


Hi Chris,

Thanks for pointing out the use of "View page source". Very helpful to know.

Do you happen to know anything about how to perform the analysis itself? I 
haven't been able to find anything confirming that the approach described in my 
original email (below) is correct. 

Thanks,

Paul


On Wed, 7/9/14, Andrews, Chris  wrote:

 Subject: RE: [R] Survival Analysis with an Historical Control
 To: "Paul Miller" , "r-help@r-project.org" 

 Received: Wednesday, July 9, 2014, 11:26 AM
 
 The code is actually
 available at the websites you provide.  Try "View page
 source" in your browser.  The most cryptic code
 isn't needed because the math functions (e.g, incomplete
 gamma function) are available in R.
 
 
 -Original Message-
 From: Paul Miller [mailto:pjmiller...@yahoo.com]
 
 Sent: Tuesday, July 08, 2014 12:00 PM
 To: r-help@r-project.org
 Subject: [R] Survival Analysis with an
 Historical Control
 
 Hello
 All,
 
 I'm trying to
 figure out how to perform a survival analysis with an
 historical control. I've spent some time looking online
 and in my boooks but haven't found much showing how to
 do this. Was wondering if there is a R package that can do
 it, or if there are resources somewhere that show the actual
 steps one takes, or if some knowledgeable person might be
 willing to share some code. 
 
 Here is a statement that describes the sort of
 analyis I'm being asked to do.
 
 A one-sample parametric test assuming an
 exponential form of survival was used to test the hypothesis
 that the treatment produces a median PFS no greater than the
 historical control PFS of 16 weeks.  A sample median PFS
 greater than 20.57 weeks would fall beyond the critical
 value associated with the null hypothesis, and would be
 considered statistically significant at alpha = .05, 1
 tailed.  
 
 My understanding
 is that the cutoff of 20.57 weeks was obtained using an
 online calculator that can be found at:
 
 http://www.swogstat.org/stat/public/one_survival.htm
 
 Thus far, I've been unable
 to determine what values were plugged into the calculator to
 get the cutoff.
 
 There's
 another calculator for a nonparamertric test that can be
 found at:
 
 http://www.swogstat.org/stat/public/one_nonparametric_survival.htm
 
 It would be nice to try doing
 this using both a parameteric and a non-parametric model.
 
 So my first question would be
 whether the approach outlined above is valid or if the
 analysis should be done some other way. If the basic idea is
 correct, is it relatively easy (for a Terry Therneau type
 genius) to implement the whole thing using R? The calculator
 is a great tool, but, if reasonable, it would be nice to be
 able to look at some code to see how the numbers actually
 get produced.
 
 Below are
 some sample survival data and code in case this proves
 helpful.
 
 Thanks,
 
 Paul
 
 ###
  Example Data: GD2 Vaccine 
 ###
 
 connection <-
 textConnection("

Re: [R] two questions - function help and 32vs64 bit sessions

[Duncan Murdoch]

On 21/07/2014 9:40 AM, Fowler, Mark wrote:

Hi Duncan,
I tried your suggestion, but no luck. The first error is no surprise, it
just confirms the address is lost. The second line suggests it worked,
but it didn't. The session is still remembering the original address.

> tools::startDynamicHelp(FALSE) # shut it down
Warning message:
In file(out, "wt") :
   cannot open file
'C:\Users\fowlerm\AppData\Local\Temp\1\RtmpK09I4B\Rhttpd26c04742884': No
such file or directory
> tools::startDynamicHelp(TRUE)  # start it up
starting httpd help server ... done


This is a different error than I thought you described.  I thought 
something had shut down the server, but it looks like something has 
deleted your temporary directory.   You might be able to tell your OS 
not to do that (if it was the OS that did), or you can put your 
temporary directory in a location that is less likely to get deleted, 
e.g. by setting the TMPDIR environment variable to somewhere else before 
you start R.  For example, I typically run with TMPDIR set to C:/temp 
when I'm debugging things.


Duncan Murdoch



-Original Message-
From: Duncan Murdoch [mailto:murdoch.dun...@gmail.com]
Sent: July 16, 2014 10:47 AM
To: Fowler, Mark; r-h...@stat.math.ethz.ch
Subject: Re: [R] two questions - function help and 32vs64 bit sessions

On 14/07/2014, 11:42 AM, Fowler, Mark wrote:
> Hello,
>
>
>
> Two unrelated questions, and neither urgent.
>
>
>
> Windows 7, R 3.0.1. Using R Console, no fancy interface.
>
>
>
> The function help ultimately becomes lost to a session kept running
> for extended periods (days). I.e. with a new session if you invoke the

> Help menu 'R functions (txt)...' it activates the html help and goes
> to the named function page. This will work fine for at least a day,
> but typically the next day invoking the help menu in this fashion will

> fail, as R looks for a temporary address it creates on your computer.
> This gets lost, possibly due to network administration activity. So
> then I save and start another session with same Rdata. Trivial enough
> but irritating. Anybody know how to restore the 'link' without ending
> and restarting the session?

I never have sessions that last that long, so I haven't tried this, but
I'd expect you could restart the help system in this way:

tools::startDynamicHelp(FALSE) # shut it down
tools::startDynamicHelp(TRUE)  # start it up

Duncan Murdoch

>
>
>
> I have a mix of 32-bit and 64-bit requirements, with 64 the default. I

> became used to starting R sessions directly from the appropriate Rdata

> workspaces. With the latest version I need to start from the generic
> icon and then load the workspace if I want 32-bit. Anybody know a way
> to make the Rdata files keep track of which bit version they work
> with, or some trick that accomplishes the same objective? The 32-bit
> requirement is usually just RODBC, and the need for it is scattered
> over lots of workspaces. Again, trivial but a nuisance. A more
> pragmatic motive is to not oblige users of applications to think about

> it. Any way to make a session switch 'bits' with a source file?
>
>
>
> Mark Fowler
> Population Ecology Division
> Bedford Inst of Oceanography
> Dept Fisheries & Oceans
> Dartmouth NS Canada
> B2Y 4A2
> Tel. (902) 426-3529
> Fax (902) 426-9710
> Email mark.fow...@dfo-mpo.gc.ca 
>
>
>
>
>
>
>[[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>



__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Semi Markov warnings ( for dummies)

[M.A. Pet]

*Hello,*

I never worked with R before my supervisor asks me to run a semiMarkov
analysis a month ago. After a long struggle, to date, the code works, but I
still get some warnings. However, because of my lack of knowledge in R I am
not possible to figure out the problems or say anything about the influence
of these warnings on my outcome. Hopefully, someone would help me with
these ( I think basic) questions.


What is the case? I want to do a semiMarkov analysis with 3 states (state
n, state s and state e). Wherefore I want to run the analysis to see
whether there is a difference in Hazard Ratio between transitions nàe and sàe.
Therefore, I’ve got data of nearly 60 persons. In excel, every worksheet
reflects a person.  To answer my research question I tried to run the
script for one person.

Unfortunately, I get the following error:

> # read data and convert to correct format (and print tail).

> data <- readdata( input <- 'E1_tijd_stress.csv' )

Error in `$<-.data.frame`(`*tmp*`, "state", value = "s") :

  replacement has 1 row, data has 0

Called from: `$<-`(`*tmp*`, "state", value = "s")



Furthermore, I used a macro in excel to merge these worksheets into one
worksheet ( all persons in one file) . However, after running the script I
saw the following warnings:

> # semi-Markov model without covariates.

> fit.semi.markov <- semiMarkov( data = data, states = states.semi.markov, 
> mtrans = mtrans.semi.markov )



Iter: 1 fn: 11224.5658  Pars:  1000.0 412.24725 261.14293
66.30789   1.67507   5.59420   0.55746   1.37715   0.59734   0.48604
0.73281   0.74205

Iter: 2 fn: 11224.5658  Pars:  1000.0 412.24425 261.14508
66.30635   1.67506   5.59420   0.55746   1.37716   0.59734   0.48604
0.73281   0.74206

solnp--> Completed in 2 iterations

Warning messages:

1: In .safefunx(tmpv, .solnp_fun, .env, ...) :

solnp-->warning: Inf detected in function call...check your function


These repeated for 7 times. Or the following warning:

> fit.semi.markov <- semiMarkov( data = data, states = states.semi.markov, 
> mtrans = mtrans.semi.markov )

Iter: 1 fn: 444562.8860 Pars:  1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0 0.47882 0.72691 0.74627

solnp--> Completed in 1 iterations

Error in svd(X) : infinite or missing values in 'x'

In addition: There were 29 warnings (use warnings() to see them)

29: In .safefunx(tmpv, .solnp_fun, .env, ...) :

solnp-->warning: NaN detected in function call...check your function


Furthermore, I want to control the analysis for some covariates.

> # semi-Markov model with covariates

> fit2 <- semiMarkov( data = data,cov= as.data.frame(data$stress.score), states 
> = states, mtrans = mtrans )

Error in mtrans[as.numeric(substring(trans.hh[i], first = 1, last = 1)),  :

  subscript out of bounds


Is someone able to tell me what these warnings exactly mean? Or someone I
can contact to discuss about this? I understand inf and NaN for example but
I am not able to figure out what to change in my data format of script (
the problem is, some of my data run without problems, some does not).
Hopefully I give enough information, otherwise let me know.

Thank you!

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] two questions - function help and 32vs64 bit sessions

[Fowler, Mark]

Hi Duncan,
I tried your suggestion, but no luck. The first error is no surprise, it
just confirms the address is lost. The second line suggests it worked,
but it didn't. The session is still remembering the original address.

> tools::startDynamicHelp(FALSE) # shut it down
Warning message:
In file(out, "wt") :
  cannot open file
'C:\Users\fowlerm\AppData\Local\Temp\1\RtmpK09I4B\Rhttpd26c04742884': No
such file or directory
> tools::startDynamicHelp(TRUE)  # start it up
starting httpd help server ... done

-Original Message-
From: Duncan Murdoch [mailto:murdoch.dun...@gmail.com] 
Sent: July 16, 2014 10:47 AM
To: Fowler, Mark; r-h...@stat.math.ethz.ch
Subject: Re: [R] two questions - function help and 32vs64 bit sessions

On 14/07/2014, 11:42 AM, Fowler, Mark wrote:
> Hello,
> 
>  
> 
> Two unrelated questions, and neither urgent.
> 
>  
> 
> Windows 7, R 3.0.1. Using R Console, no fancy interface.
> 
>  
> 
> The function help ultimately becomes lost to a session kept running 
> for extended periods (days). I.e. with a new session if you invoke the

> Help menu 'R functions (txt)...' it activates the html help and goes 
> to the named function page. This will work fine for at least a day, 
> but typically the next day invoking the help menu in this fashion will

> fail, as R looks for a temporary address it creates on your computer. 
> This gets lost, possibly due to network administration activity. So 
> then I save and start another session with same Rdata. Trivial enough 
> but irritating. Anybody know how to restore the 'link' without ending 
> and restarting the session?

I never have sessions that last that long, so I haven't tried this, but
I'd expect you could restart the help system in this way:

tools::startDynamicHelp(FALSE) # shut it down
tools::startDynamicHelp(TRUE)  # start it up

Duncan Murdoch

> 
>  
> 
> I have a mix of 32-bit and 64-bit requirements, with 64 the default. I

> became used to starting R sessions directly from the appropriate Rdata

> workspaces. With the latest version I need to start from the generic 
> icon and then load the workspace if I want 32-bit. Anybody know a way 
> to make the Rdata files keep track of which bit version they work 
> with, or some trick that accomplishes the same objective? The 32-bit 
> requirement is usually just RODBC, and the need for it is scattered 
> over lots of workspaces. Again, trivial but a nuisance. A more 
> pragmatic motive is to not oblige users of applications to think about

> it. Any way to make a session switch 'bits' with a source file?
> 
>  
> 
> Mark Fowler
> Population Ecology Division
> Bedford Inst of Oceanography
> Dept Fisheries & Oceans
> Dartmouth NS Canada
> B2Y 4A2
> Tel. (902) 426-3529
> Fax (902) 426-9710
> Email mark.fow...@dfo-mpo.gc.ca 
> 
> 
> 
>  
> 
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
> 

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] Include plotting symbals pch 16 and 17 into captions / text in pdf graph

[Rainer M Krug]

Hi

I want to include the plotting symbols pch=16 and pch=17 in text in
the graph used as labels.

This works:

--8<---cut here---start->8---
plot(0:1,0:1, type="n")
points(0.2, 1, pch=16)
mtext(
text = "pch=16 (\U25CF)",
side = 3,
at   = 0.2,
line = 1,
)
points(0.8, 1, pch=17)
mtext(
text = "pch=17 (\U25B2)",
side = 3,
at   = 0.8,
line = 1
)
--8<---cut here---end--->8---

But I would need this as a pdf and in Calibri font, which does not
include the UTF symbol \U25B2 (the triangle). I got that far:

--8<---cut here---start->8---

cairo_pdf("utftext.pdf", family="Calibri")
plot(0:1,0:1, type="n")
points(0.2, 1, pch=16)
mtext(
text = "pch=16 (\U25CF)",
side = 3,
at   = 0.2,
line = 1,
)
points(0.8, 1, pch=17)
mtext(
text = "pch=17 (\U25B2)",
side = 3,
at   = 0.8,
line = 1
)
dev.off()
--8<---cut here---end--->8---

which results in the attached pdf (hope it comes through). If not: the
utf symbol for the upward error is displayed as an empty square with a
question mark in it).

Is there a way that I can show the plotting symbol pch=17 in the caption
when using this font?

Thanks,

Rainer

-- 
Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology, 
UCT), Dipl. Phys. (Germany)

Centre of Excellence for Invasion Biology
Stellenbosch University
South Africa

Tel :   +33 - (0)9 53 10 27 44
Cell:   +33 - (0)6 85 62 59 98
Fax :   +33 - (0)9 58 10 27 44

Fax (D):+49 - (0)3 21 21 25 22 44

email:  rai...@krugs.de

Skype:  RMkrug

PGP: 0x0F52F982


utftext.pdf
Description: Adobe PDF document


pgpfeJIhfMA2r.pgp
Description: PGP signature
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] packages across different versions of R

[Prof Brian Ripley]

On 21/07/2014 12:24, Charles Thuo wrote:

  I have just installed R 3.1.1 in a machine where R 3.0.1 is already
installed.  Is it possible to use packages in the 3.0.1 on the 3.1.1.
version as the same are in a single workstation.


Perhaps, perhaps not.  It depends in part on your platform which you 
have not told us (see the posting guide).  It is definitely safer to 
reinstall them.  For Window users see 
http://cran.r-project.org/bin/windows/base/rw-FAQ.html#What_0027s-the-best-way-to-upgrade_003f 
: similar ideas work for other platforms.


The place where using packages for 3.0.1 is least likely to work is OS 
X, but there installing a new binary version of R by default removes 
earlier installations.





Charles

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.




--
Brian D. Ripley,  rip...@stats.ox.ac.uk
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@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] packages across different versions of R

[Charles Thuo]

 I have just installed R 3.1.1 in a machine where R 3.0.1 is already
installed.  Is it possible to use packages in the 3.0.1 on the 3.1.1.
version as the same are in a single workstation.


Charles

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

Re: [R] plotly

[Sarah Goslee]

Hi,

On Monday, July 21, 2014, Shane Carey  wrote:

> Hey,
>
> What version of R is required to use the plotly library?
>
> I have R version 3.0.1 and it will not allow me to install the devtools
> package or the ploty package.
>
> I have googled and searched to see what version of R I should be running
> but could not find anything.


It's always a good idea to upgrade to the current release of R before
asking questions like that (3.1.1). But in general, packages on CRAN
clearly state what version of R is needed, as in

http://cran.r-project.org/web/packages/devtools/index.html
devtools: Tools to make developing R code easier

Collection of package development tools

Version:1.5Depends:R (≥ 3.0.2)

For packages not on CRAN, like plotly, you may need to download the package
and check the DESCRIPTION file.

Sarah


-- 
Sarah Goslee
http://www.stringpage.com
http://www.sarahgoslee.com
http://www.functionaldiversity.org

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[R] plotly

[Shane Carey]

Hey,

What version of R is required to use the plotly library?

I have R version 3.0.1 and it will not allow me to install the devtools
package or the ploty package.

I have googled and searched to see what version of R I should be running
but could not find anything.

Thanks

-- 
Shane

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.