[R] Plot prediction for Heckman 2-step (sampleSelection) models after predict.selection

[Syed Rashid Munir]

Hi everyone,

I'm fairly new to R and I've run into a problem that I can't find a 
solution for.

I'm trying to plot the predicted values of a dependent variable based on 
specific values of one independent variable after running a Heckit 
(Heckman 2-step) model from the sampleSelection package.

The general model looks like this:

|Mod1 <- Heckit(selection = binaryindicator ~ x1 + z, outcome = y ~ x1, 
df, method) |

where, the data structure is as follows:

|z x1 y binaryindicator 0.5 4 300 1 0.1 8 400 1 0.2 10 500 1 0.2 18 NA 0 
0.4 20 50 1 0.3 30 1000 1 |

After pouring over documentation for sampleSelection, I can only extract 
predictions from the predict.selection function and get the 
corresponding residuals. However, this method doesn't let me specify 
predictions for one variable (say x1) and resulting predictions are from 
the entire model.

Therefore, I am looking for a way (canned function or manual) to plot 
the prediction on y from x1 with confidence intervals. All attempts to 
do so (using plot, cplot, ggpredict, effects_plot) have failed. All of 
the above functions work with linear models, but not a selection model.

I'd ideally want to also specify which equation (outcome or selection) I 
want the prediction from, but that has also eluded me, as I can't get 
any predictions to begin with. Is there a way to achieve this? In Stata, 
this can be done with the margins command, but the equivalent margins 
package in R doesn't work for sampleSelection models (for me so far). 
Any and all help would be much appreciated!


Regards,
Rashid
Poli Sci Student, SUNY Binghamton

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] converting zipcodes to latitude/longitude

[Jim Lemon]

Hi Nicola,
Getting the blank rows will be a bit more difficult and I don't see
why they should be in the final data frame, so:

townzip<-read.table(text="waltham, Massachusetts 02451
Columbia, SC 29209

Wheat Ridge , Colorado 80033
Charlottesville, Virginia 22902
Fairbanks, AK 99709
Montpelier, VT 05602
Dobbs Ferry, New York 10522

Henderson , Kentucky 42420",
sep="\t",stringsAsFactors=FALSE)
zip_split<-function(x) {
 commasplit<-unlist(strsplit(x,","))
 state<-trimws(gsub("[[:digit:]]","",commasplit[2]))
 zip<-trimws(gsub("[[:alpha:]]","",commasplit[2]))
 return(c(commasplit[1],state,zip))
}
townzipsplit<-as.data.frame(t(sapply(townzip$V1,zip_split)))
rownames(townzipsplit)<-NULL
names(townzipsplit)<-c("town","state","zip")
townzipsplit$latlon<-NA
# I don't know the name of the zipcode column in the "zipcode" data frame
newzipdf<-merge(townzipsplit,zipcodedf,by.x="zip",by.y="zip")

Jim

On Tue, May 14, 2019 at 5:57 AM Nicola Ruggiero
 wrote:
>
> Hello everyone,
>
> I've downloaded Jeffrey Breen's R package "zipcode," which has the
> latitude and longitude for all of the US zip codes. So, this is a
> data.frame with 43,191 observations. That's one data frame in my
> environment.
>
> Then, I have another data.frame with over 100,000 observations that
> look like this:
>
> waltham, Massachusetts 02451
> Columbia, SC 29209
>
> Wheat Ridge , Colorado 80033
> Charlottesville, Virginia 22902
> Fairbanks, AK 99709
> Montpelier, VT 05602
> Dobbs Ferry, New York 10522
>
> Henderson , Kentucky 42420
>
> The spaces represent absences in the column. Regardless,
> I need to figure out how to write a code that would, presumably, match
> the zipcodes and produce another column to the data frame with the
> latitude and longitude. So, for example, the code would recognize
> 02451 above, and, in the the column next to it, the code would write
> 42.3765° N, 71.2356° W in the column next to it, since that's the
> latitude and longitude for Waltham, Massachusetts.
>
> Any idea of how to begin a code that would perform such an operation?
>
> Again, I have a data.frame with the zipcodes linked to the the
> latitudes and longitudes, on the one hand, and another data.frame with
> only zipcodes (and some holes). I need to produce the corresponding
> latitude/longitudes in the latter data.frame.
>
> Nicola
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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 -- To UNSUBSCRIBE and more, see
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] install pcre

[Jeff Newmiller]

Looks like you don't have all of the dev tools installed on your unnamed OS. I 
highly recommend Googling for how to install R in your OS.

Please read the Posting Guide... Discussions about compiling R should be in 
R-devel. Most people install binaries which would be relevant here.

On May 13, 2019 8:31:22 PM MDT, yueli  wrote:
>Hello,
>
>I am trying to install R.
>
>Thanks in advance for any help!
>
>Yue
>
>
>checking for pcre.h... yes
>checking pcre/pcre.h usability... no
>checking pcre/pcre.h presence... no
>checking for pcre/pcre.h... no
>checking if PCRE version >= 8.20, < 10.0 and has UTF-8 support... no
>checking whether PCRE support suffices... configure: error: pcre >=
>8.20 library and headers are required
>
>   [[alternative HTML version deleted]]
>
>__
>R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>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.

-- 
Sent from my phone. Please excuse my brevity.

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] install pcre

[yueli]

Hello,

I am trying to install R.

Thanks in advance for any help!

Yue


checking for pcre.h... yes
checking pcre/pcre.h usability... no
checking pcre/pcre.h presence... no
checking for pcre/pcre.h... no
checking if PCRE version >= 8.20, < 10.0 and has UTF-8 support... no
checking whether PCRE support suffices... configure: error: pcre >= 8.20 
library and headers are required

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Trying to understand the magic of lm (Still trying)

[Sorkin, John]

Terry ,


Thank you. Many years ago I took a course you taught in which you explained how 
to conduct survival analyses using S. The course was very useful, as was the 
email you sent me today. If you find a place where you can store your lecture 
notes, please send me the URL.


I believe that there is a great need for someone to explain not just how to 
write a package, but generally how to write a function that checks the 
parameters passed to it, that uses the parameters in a manner that allows the 
output of the function to produce output that informs the function users of the 
call to the function, etc. While these steps are needed when one writes a 
package, they should be taught as a matter of good coding practice when anyone 
writes a function that will be used more than once. Many years ago when I was a 
mainframe system programmer, it was de rigueur that one learned (and used) 
certain standards about saving registers at the beginning of a function and 
restoring them at the end of the function. The same should be true for all R 
functions; certain standardized, well described steps should be considered a 
part of writing any function. The problem, at least from my perspective, is 
that there is no commonly recognized document that explains the steps clearly.


Thank you,

John


John David Sorkin M.D., Ph.D.
Professor of Medicine
Chief, Biostatistics and Informatics
University of Maryland School of Medicine Division of Gerontology and Geriatric 
Medicine
Baltimore VA Medical Center
10 North Greene Street
GRECC (BT/18/GR)
Baltimore, MD 21201-1524
(Phone) 410-605-7119
(Fax) 410-605-7913 (Please call phone number above prior to faxing)




From: Therneau, Terry M., Ph.D. 
Sent: Monday, May 13, 2019 9:29 AM
To: r-help@r-project.org; Sorkin, John
Subject: Re: [R] Trying to understand the magic of lm (Still trying)

John,

 The text below is cut out of a "how to write a package" course I gave at the R 
conference in Vanderbilt.   I need to find a home for the course notes, because 
it had a lot of tidbits that are not well explained in the R documentation.
Terry T.



Model frames:
One of the first tasks of any modeling routine is to construct a special data 
frame containing the covariates that will be used, via a call to the 
model.frame function. The code to do this is found in many routines, and can be 
a little opaque on first view. The obvious code would be
\begin{verbatim}
coxph <- function(formula, data, weights, subset, na.action,
init, control, ties= c("efron", "breslow", "exact"),
singular.ok =TRUE, robust=FALSE,
model=FALSE, x=FALSE, y=TRUE,  tt, method=ties, ...) {

 mf <- model.frame(formula, data, subset, weights, na.action)
\end{verbatim}
since those are the coxph arguments that are passed forward to the model.frame 
routine.  However, this simple approach will fail with a ``not found'' error 
message if any of the data, subset, weights, etc. arguments are missing. 
Programs have to take the slightly more complicated approach of constructing a 
call.
\begin{verbatim}
Call <- match.call()
indx <- match(c("formula", "data", "weights", "subset", "na.action"),
  names(Call), nomatch=0)
if (indx[1] ==0) stop("A formula argument is required")
temp <- Call[c(1,indx)]  # only keep the arguments we wanted
temp[[1]] <- as.name('model.frame')  # change the function called
mf <- eval(temp, parent.frame())

Y <- model.response(mf)
etc.
\end{verbatim}

We start with a copy of the call to the program, which we want to save anyway 
as documentation in the output object. Then subscripting is used to extract 
only the portions of the call that we want, saving the result in a temporary. 
This is based on the fact that a call object can be viewed as a list whose 
first element is the name of the function to call, followed by the arguments to 
the call. Note the use of \code{nomatch=0}; if any arguments on the list are 
missing they will then be missing in \code{temp}, without generating an error 
message. The \mycode{temp} variable will contain a object of type ``call'', 
which is an unevaluated call to a routine.  Finally, the name of the function 
to be called is changed from ``coxph'' to ``model.frame'' and the call is 
evaluated.  In many of the core routines the result is stored in a variable 
``m''.  This is a horribly short and non-descriptive name. (The above used mf 
which isn't a much better.)  Many routines also use ``m'' for the temporar
 y variable leading to \code{m <- eval(m, parent.frame())}, but I think that is 
unnecessarily confusing.

The list of names in the match call will include all arguments that should be 
evaluated within context of the named dataframe. This can include more than the 
list above, the survfit routine for instance has an optional argument ``id'' 
that names an identifying variable (several rows of the data may represent a 
single subject), and this is included along with ``formula'' 

Re: [R-es] Como reordenar datos para analisis multiples correspondencias (MCA)

[Juan Abasolo]

Gracias, Javier;
Anduve buscando entre tus dos propuestas para mi encarables (plyr y
reshape2) y presupongo que el problema es eso del ingenio que comentás, que
hay que tener y no tengo.

Los paquetes me vinieron igualmente re-bien, porque me di cuenta que no
tengo por qué sufrir haciéndolo a mano cada vez que encaro algunas de esas
situaciones.

Pero lo de ordenar los datos para el MCA... el ingenio estara en alguna
inlampara de casa, porque a mí no se me enciende la bombilla. Parece
prudente, a falta de más, idea purgar las respuestas dobles según algún
criterio, que siempre es deformar la realidad un poquito más.



Hau idatzi du Javier Marcuzzi (javier.ruben.marcu...@gmail.com)
erabiltzaileak (2019 mai. 13, al. (14:47)):

> Estimado Juan Abasolo
>
> Para esas actividades se necesita algo de ingenio, no hay una sola
> solución, para esto existen librerías como plyr, reshape2, entre otras, o
> simplemente a mano, en lo personal mis primeros pasos para ordenar fue
> sqldf porque conocía sql más que R, sin embargo hoy las alternativas al
> respecto son muy amplias.
>
> Yo pienso en como debo colocar los datos en el modelo, luego en como tengo
> los datos originales, y finalmente busco la forma de acomodarlos. En este
> paso posiblemente es donde R tiene la mayor cantidad de alternativas e
> insumo de tiempo razonando puesto que tiene algo de "artesanal en R".
>
> Javier Rubén Marcuzzi
>
> El dom., 12 may. 2019 a las 18:53, Juan Abasolo ()
> escribió:
>
>> Necesito luz para ordenar unos datos... en realidad, para ordenar muchos
>> muchoas veces. El problema hoy:
>>
>> Tengo una base de datos que incluye respuestas dobles en algunas variables
>> en algunos individuos, tabla didáctica:
>>
>> idioma   alergia  color
>> individuo1  en,es 0   amarillo
>> individuo2  es,en huevo   limon
>> individuo3  es,fr,en  pescado, huevo  salmon
>>
>> Necesitaría ordenarlos de otra manera (creo), para poder hacer un análisis
>> de multiples concordancias, porque tal y como lo tengo puesto, por
>> ejemplo,
>> en idioma no hay concordancia, aunque es evidente que entre todos podrían
>> comunicarse entre castellano o inglés, por ejemplo o que son equivalentes,
>> también en idioma, los sujetos 1 y 2.
>>
>> Los datos con los que tengo que trabajar, de momento, son equivalentes a
>> los del ejemplo, lease categóricos e incluyen posibilidad de respuestas
>> múltiples.
>>
>> Supongo que le resultará obvio a alguno cómo resolverlo o en qué está mal
>> el orden así de los datos... pero yo nopuedo darme cuenta.
>>
>> a) Alguna pista?
>> b) Recomendación de lectura para abrir un poco la mente? (tengo más de una
>> de este estilo)
>>
>> Gracias
>>
>>
>>
>> --
>> Juan Abasolo
>>
>> Hizkuntzaren eta Literaturaren Didaktika Saila | EUDIA ikerketa taldea
>> Bilboko Hezkuntza Fakultatea
>> Euskal Herriko Unibertsitatea
>> UPV/EHU
>>
>> Sarriena auzoa z/g 48940 - Leioa (Bizkaia)
>>
>> T: (+34) 94 601 7567
>> Telegram: @JuanAbasolo
>> Skype: abasolo72
>>
>> Tutoretza ordutegia 
>>
>> [[alternative HTML version deleted]]
>>
>> ___
>> R-help-es mailing list
>> R-help-es@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-help-es
>>
>

-- 
Juan Abasolo

Hizkuntzaren eta Literaturaren Didaktika Saila | EUDIA ikerketa taldea
Bilboko Hezkuntza Fakultatea
Euskal Herriko Unibertsitatea
UPV/EHU

Sarriena auzoa z/g 48940 - Leioa (Bizkaia)

T: (+34) 94 601 7567
Telegram: @JuanAbasolo
Skype: abasolo72

Tutoretza ordutegia 

[[alternative HTML version deleted]]

___
R-help-es mailing list
R-help-es@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es

[R] converting zipcodes to latitude/longitude

[Nicola Ruggiero]

Hello everyone,

I've downloaded Jeffrey Breen's R package "zipcode," which has the
latitude and longitude for all of the US zip codes. So, this is a
data.frame with 43,191 observations. That's one data frame in my
environment.

Then, I have another data.frame with over 100,000 observations that
look like this:

waltham, Massachusetts 02451
Columbia, SC 29209

Wheat Ridge , Colorado 80033
Charlottesville, Virginia 22902
Fairbanks, AK 99709
Montpelier, VT 05602
Dobbs Ferry, New York 10522

Henderson , Kentucky 42420

The spaces represent absences in the column. Regardless,
I need to figure out how to write a code that would, presumably, match
the zipcodes and produce another column to the data frame with the
latitude and longitude. So, for example, the code would recognize
02451 above, and, in the the column next to it, the code would write
42.3765° N, 71.2356° W in the column next to it, since that's the
latitude and longitude for Waltham, Massachusetts.

Any idea of how to begin a code that would perform such an operation?

Again, I have a data.frame with the zipcodes linked to the the
latitudes and longitudes, on the one hand, and another data.frame with
only zipcodes (and some holes). I need to produce the corresponding
latitude/longitudes in the latter data.frame.

Nicola

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Error message when adding drift for Arima model

[Michael Howell]

Yes that is a little off topic but I will look into it more. Thank you very
much for your help.

Michael

On Mon, May 13, 2019 at 11:33 AM Rui Barradas  wrote:

> Hello,
>
> Sorry for the late reply.
> Inline.
>
> Às 17:54 de 10/05/19, Michael Howell escreveu:
> > Rui,
> > I'm still new to ARIMA forecasting but I examined the PACF and saw
> > significant correlation at lag 2.
>
> You saw a PACF with a significant correlation at lag 2 but not at lag 1.
> When this happens, it many times means that you shouldn't consider the
> lag 2. In fact, it might mean that the process is nonlinear.
> And the ACF shows an insignificant lag 1.
>
> Try
>
> ords <- list(c(1, 0, 0), c(2, 0, 0), c(0, 0, 1))
> fit_list <- lapply(ords, function(o)
>Arima(tsdata, order = o, include.drift = TRUE))
> sapply(fit_list, AIC)
> sapply(fit_list, BIC)
>
>
> Which gives the minimum AIC? And BIC?
> These are not perfect and automated model selection can have problems,
> but it's not unreasonable to compare them.
>
> I believe this is off-topic for R-Help, since it's a question about
> statistics and nonlinear time series is a really, really broad field to
> be discussed here. Try to find local help on this.
>
> Hope this helps,
>
> Rui Barradas
>
> The ACF showed a more gradual decline
> > which seemed to indicate it was Autoregressive. That should mean it's a
> > AR(2) process right?
> >
> > image.png
> > **//___^
> > Regards,
> > Michael Howell
> >
> >
> > On Fri, May 10, 2019 at 12:51 AM Rui Barradas  > > wrote:
> >
> > Why not
> >
> > Arima(tsdata, c(0, 0, 1), include.drift = TRUE)
> >
> > ?
> >
> > Why do you say it should be an AR(2) model?
> >
> > Hope this helps,
> >
> > Rui Barradas
> >
> > Às 06:43 de 10/05/19, Rui Barradas escreveu:
> >  > Hello,
> >  >
> >  > This is just a typo, in R logical values ("true) are not character
> >  > strings. You must pass FALSE (the default, can be omited) or TRUE.
> >  >
> >  > fitdata  <-  Arima(tsdata, c(2, 0, 0), include.drift = TRUE)
> >  >
> >  >
> >  >  From the help page ?logical
> >  >
> >  > Details
> >  >
> >  > TRUE and FALSE are reserved words denoting logical constants in
> > the R
> >  > language, whereas T and F are global variables whose initial
> > values set
> >  > to these. All four are logical(1) vectors.
> >  >
> >  > Hope this helps,
> >  >
> >  > Rui Barradas
> >  >
> >  > Às 00:26 de 10/05/19, Bert Gunter escreveu:
> >  >> In future, always cc the list (unless it's personal,which this
> > isn't). I
> >  >> have done so here. As I am largely ignorant on the subject
> > matter, others
> >  >> will have to help, which is why you should cc the list.
> >  >>
> >  >> Cheers,
> >  >> Bert Gunter
> >  >>
> >  >> "The trouble with having an open mind is that people keep coming
> > along
> >  >> and
> >  >> sticking things into it."
> >  >> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip
> )
> >  >>
> >  >>
> >  >> On Thu, May 9, 2019 at 3:49 PM Michael Howell
> > mailto:mchowe...@gmail.com>>
> >  >> wrote:
> >  >>
> >  >>> I apologize for that. The Arima() function that I'm trying to
> > use comes
> >  >>> from the forecast package. I created a time series object using
> > the
> >  >>> above
> >  >>> 24 observations. The initial model I created doesn't seem to
> > perform so
> >  >>> well so I thought a drift term might fit the data better. I
> > used the
> >  >>> following code to create the time series object:
> >  >>>
> >  >>> tsdata<- ts(data, start = c(1,1), end = c(24,1), frequency = 1)
> >  >>>
> >  >>>
> >  >>> Where* data* is the dataframe that contains the initial 24
> > observations.
> >  >>> I then used the following code to try to create the model:
> >  >>>
> >  >>> fitdata  <-  Arima(tsdata,c(2,0,0),include.drift="true")
> >  
> >  >>>
> >  >>> After doing this I obtained the following error message:
> >  >>>
> >  >>> Error in (order[2] + seasonal$order[2]) > 1 & include.drift:
> > operations
> >   are possible only for numeric, logical or complex types
> >   Traceback:
> >  
> >   1. Arima(tsdata, c(2, 0, 0), include.drift = "true")
> >  >>>
> >  >>>
> >  >>>   I hope this is more clear.
> >  >>>
> >  >>> On Thu, May 9, 2019 at 4:39 PM Bert Gunter
> > mailto:bgunter.4...@gmail.com>>
> >  >>> wrote:
> >  >>>
> >   Please start by reading and following the posting guide linked
> > at the
> >   bottom of this email. In particular:
> >  
> >   1) Post in **plain text** on this plain text list so we don't
> > get the
> >   mangled html of your post.
> >  
> >   2) Tell us what package Arima() is 

Re: [R] Error message when adding drift for Arima model

[Rui Barradas]

Hello,

Sorry for the late reply.
Inline.

Às 17:54 de 10/05/19, Michael Howell escreveu:

Rui,
I'm still new to ARIMA forecasting but I examined the PACF and saw 
significant correlation at lag 2. 


You saw a PACF with a significant correlation at lag 2 but not at lag 1. 
When this happens, it many times means that you shouldn't consider the 
lag 2. In fact, it might mean that the process is nonlinear.

And the ACF shows an insignificant lag 1.

Try

ords <- list(c(1, 0, 0), c(2, 0, 0), c(0, 0, 1))
fit_list <- lapply(ords, function(o)
  Arima(tsdata, order = o, include.drift = TRUE))
sapply(fit_list, AIC)
sapply(fit_list, BIC)


Which gives the minimum AIC? And BIC?
These are not perfect and automated model selection can have problems, 
but it's not unreasonable to compare them.


I believe this is off-topic for R-Help, since it's a question about 
statistics and nonlinear time series is a really, really broad field to 
be discussed here. Try to find local help on this.


Hope this helps,

Rui Barradas

The ACF showed a more gradual decline
which seemed to indicate it was Autoregressive. That should mean it's a 
AR(2) process right?


image.png
**//___^
Regards,
Michael Howell


On Fri, May 10, 2019 at 12:51 AM Rui Barradas > wrote:


Why not

Arima(tsdata, c(0, 0, 1), include.drift = TRUE)

?

Why do you say it should be an AR(2) model?

Hope this helps,

Rui Barradas

Às 06:43 de 10/05/19, Rui Barradas escreveu:
 > Hello,
 >
 > This is just a typo, in R logical values ("true) are not character
 > strings. You must pass FALSE (the default, can be omited) or TRUE.
 >
 > fitdata  <-  Arima(tsdata, c(2, 0, 0), include.drift = TRUE)
 >
 >
 >  From the help page ?logical
 >
 > Details
 >
 > TRUE and FALSE are reserved words denoting logical constants in
the R
 > language, whereas T and F are global variables whose initial
values set
 > to these. All four are logical(1) vectors.
 >
 > Hope this helps,
 >
 > Rui Barradas
 >
 > Às 00:26 de 10/05/19, Bert Gunter escreveu:
 >> In future, always cc the list (unless it's personal,which this
isn't). I
 >> have done so here. As I am largely ignorant on the subject
matter, others
 >> will have to help, which is why you should cc the list.
 >>
 >> Cheers,
 >> Bert Gunter
 >>
 >> "The trouble with having an open mind is that people keep coming
along
 >> and
 >> sticking things into it."
 >> -- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
 >>
 >>
 >> On Thu, May 9, 2019 at 3:49 PM Michael Howell
mailto:mchowe...@gmail.com>>
 >> wrote:
 >>
 >>> I apologize for that. The Arima() function that I'm trying to
use comes
 >>> from the forecast package. I created a time series object using
the
 >>> above
 >>> 24 observations. The initial model I created doesn't seem to
perform so
 >>> well so I thought a drift term might fit the data better. I
used the
 >>> following code to create the time series object:
 >>>
 >>> tsdata<- ts(data, start = c(1,1), end = c(24,1), frequency = 1)
 >>>
 >>>
 >>> Where* data* is the dataframe that contains the initial 24
observations.
 >>> I then used the following code to try to create the model:
 >>>
 >>> fitdata  <-  Arima(tsdata,c(2,0,0),include.drift="true")
 
 >>>
 >>> After doing this I obtained the following error message:
 >>>
 >>> Error in (order[2] + seasonal$order[2]) > 1 & include.drift:
operations
  are possible only for numeric, logical or complex types
  Traceback:
 
  1. Arima(tsdata, c(2, 0, 0), include.drift = "true")
 >>>
 >>>
 >>>   I hope this is more clear.
 >>>
 >>> On Thu, May 9, 2019 at 4:39 PM Bert Gunter
mailto:bgunter.4...@gmail.com>>
 >>> wrote:
 >>>
  Please start by reading and following the posting guide linked
at the
  bottom of this email. In particular:
 
  1) Post in **plain text** on this plain text list so we don't
get the
  mangled html of your post.
 
  2) Tell us what package Arima() is in.
 
  Cheers,
  Bert Gunter
 
 
 
 
  On Thu, May 9, 2019 at 2:27 PM Michael Howell
mailto:mchowe...@gmail.com>>
  wrote:
 
 > Hello everyone,
 > So this is my first post to this list, I'm trying to fit an
Arima
 > (2,0,0)
 > model and I think a drift term would help but I'm getting an
error
 > term
 > when I'm trying to include it. Here is my data:
 >
 > -6.732172338
 > -2.868884273
 > -5.371585089
 > -6.512740463
 > -4.171062657
 > -5.738499071
 > 

Re: [R] linear model contrast in R

[Richard M. Heiberger]

I think you might be looking for
?contrasts
to form the contrast matrix.

Rich

On Mon, May 13, 2019 at 7:31 AM Witold E Wolski  wrote:
>
> I am looking for a function to compute contrasts with a interface
> similar to that of
>
> lmerTest::contest
> multcomp::glht
>
> i.e. taking the model and a contrast vector or matrix as an argument,
> but for linear models, and without the multiple testing adjusted made
> by multcomp::glht.
>
> Thank you
>
>
> --
> Witold Eryk Wolski
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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 -- To UNSUBSCRIBE and more, see
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] Trying to understand the magic of lm (Still trying)

[Therneau, Terry M., Ph.D. via R-help]

John,

  The text below is cut out of a "how to write a package" course I gave at the 
R 
conference in Vanderbilt.   I need to find a home for the course notes, because 
it had a 
lot of tidbits that are not well explained in the R documentation.
Terry T.



Model frames:
One of the first tasks of any modeling routine is to construct a special data 
frame 
containing the covariates that will be used, via a call to the model.frame 
function. The 
code to do this is found in many routines, and can be a little opaque on first 
view. The 
obvious code would be
\begin{verbatim}
coxph <- function(formula, data, weights, subset, na.action,
     init, control, ties= c("efron", "breslow", "exact"),
     singular.ok =TRUE, robust=FALSE,
     model=FALSE, x=FALSE, y=TRUE,  tt, method=ties, ...) {

  mf <- model.frame(formula, data, subset, weights, na.action)
\end{verbatim}
since those are the coxph arguments that are passed forward to the model.frame 
routine.  
However, this simple approach will fail with a ``not found'' error message if 
any of the 
data, subset, weights, etc. arguments are missing. Programs have to take the 
slightly more 
complicated approach of constructing a call.
\begin{verbatim}
Call <- match.call()
indx <- match(c("formula", "data", "weights", "subset", "na.action"),
   names(Call), nomatch=0)
if (indx[1] ==0) stop("A formula argument is required")
temp <- Call[c(1,indx)]  # only keep the arguments we wanted
temp[[1]] <- as.name('model.frame')  # change the function called
mf <- eval(temp, parent.frame())

Y <- model.response(mf)
etc.
\end{verbatim}

We start with a copy of the call to the program, which we want to save anyway 
as 
documentation in the output object. Then subscripting is used to extract only 
the portions 
of the call that we want, saving the result in a temporary. This is based on 
the fact that 
a call object can be viewed as a list whose first element is the name of the 
function to 
call, followed by the arguments to the call. Note the use of \code{nomatch=0}; 
if any 
arguments on the list are missing they will then be missing in \code{temp}, 
without 
generating an error message. The \mycode{temp} variable will contain a object 
of type 
``call'', which is an unevaluated call to a routine.  Finally, the name of the 
function to 
be called is changed from ``coxph'' to ``model.frame'' and the call is 
evaluated.  In many 
of the core routines the result is stored in a variable ``m''.  This is a 
horribly short 
and non-descriptive name. (The above used mf which isn't a much better.)  Many 
routines 
also use ``m'' for the temporary variable leading to \code{m <- eval(m, 
parent.frame())}, 
but I think that is unnecessarily confusing.

The list of names in the match call will include all arguments that should be 
evaluated 
within context of the named dataframe. This can include more than the list 
above, the 
survfit routine for instance has an optional argument ``id'' that names an 
identifying 
variable (several rows of the data may represent a single subject), and this is 
included 
along with ``formula'' etc in the list of choices in the match function.  The 
order of 
names in the list makes no difference.  The id is later retrieved with 
\code{model.extract(m, 'id')}, which will be NULL if the argument was not 
supplied. At the 
time that coxph was written I had not caught on to this fact and thought that 
all 
variables that came from a data frame had to be represented in the formula 
somehow, thus 
the use of \code{cluster(id)} as part of the formula, in order to denote a 
grouping variable.

On 5/11/19 5:00 AM, r-help-requ...@r-project.org wrote:
> A number of people have helped me in my mission to understand how lm (and 
> other fucntions) are able to pass a dataframe and then refer to a specific 
> column in the dataframe. I thank everyone who has responded. I now know a bit 
> about deparse(substitute(xx)), but I still don't fully understand how it 
> works. The program below attempts to print a column of a dataframe from a 
> function whose parameters include the dataframe (df) and the column requested 
> (col). The program works fine until the last print statement were I receive 
> an error,  Error in `[.data.frame`(df, , col) : object 'y' not found . I hope 
> someone can explain to me (1) why my code does not work, and (2) what I can 
> do to fix it.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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-es] Como reordenar datos para analisis multiples correspondencias (MCA)

[Javier Marcuzzi]

Estimado Juan Abasolo

Para esas actividades se necesita algo de ingenio, no hay una sola
solución, para esto existen librerías como plyr, reshape2, entre otras, o
simplemente a mano, en lo personal mis primeros pasos para ordenar fue
sqldf porque conocía sql más que R, sin embargo hoy las alternativas al
respecto son muy amplias.

Yo pienso en como debo colocar los datos en el modelo, luego en como tengo
los datos originales, y finalmente busco la forma de acomodarlos. En este
paso posiblemente es donde R tiene la mayor cantidad de alternativas e
insumo de tiempo razonando puesto que tiene algo de "artesanal en R".

Javier Rubén Marcuzzi

El dom., 12 may. 2019 a las 18:53, Juan Abasolo ()
escribió:

> Necesito luz para ordenar unos datos... en realidad, para ordenar muchos
> muchoas veces. El problema hoy:
>
> Tengo una base de datos que incluye respuestas dobles en algunas variables
> en algunos individuos, tabla didáctica:
>
> idioma   alergia  color
> individuo1  en,es 0   amarillo
> individuo2  es,en huevo   limon
> individuo3  es,fr,en  pescado, huevo  salmon
>
> Necesitaría ordenarlos de otra manera (creo), para poder hacer un análisis
> de multiples concordancias, porque tal y como lo tengo puesto, por ejemplo,
> en idioma no hay concordancia, aunque es evidente que entre todos podrían
> comunicarse entre castellano o inglés, por ejemplo o que son equivalentes,
> también en idioma, los sujetos 1 y 2.
>
> Los datos con los que tengo que trabajar, de momento, son equivalentes a
> los del ejemplo, lease categóricos e incluyen posibilidad de respuestas
> múltiples.
>
> Supongo que le resultará obvio a alguno cómo resolverlo o en qué está mal
> el orden así de los datos... pero yo nopuedo darme cuenta.
>
> a) Alguna pista?
> b) Recomendación de lectura para abrir un poco la mente? (tengo más de una
> de este estilo)
>
> Gracias
>
>
>
> --
> Juan Abasolo
>
> Hizkuntzaren eta Literaturaren Didaktika Saila | EUDIA ikerketa taldea
> Bilboko Hezkuntza Fakultatea
> Euskal Herriko Unibertsitatea
> UPV/EHU
>
> Sarriena auzoa z/g 48940 - Leioa (Bizkaia)
>
> T: (+34) 94 601 7567
> Telegram: @JuanAbasolo
> Skype: abasolo72
>
> Tutoretza ordutegia 
>
> [[alternative HTML version deleted]]
>
> ___
> R-help-es mailing list
> R-help-es@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-help-es
>

[[alternative HTML version deleted]]

___
R-help-es mailing list
R-help-es@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-help-es

[R] linear model contrast in R

[Witold E Wolski]

I am looking for a function to compute contrasts with a interface
similar to that of

lmerTest::contest
multcomp::glht

i.e. taking the model and a contrast vector or matrix as an argument,
but for linear models, and without the multiple testing adjusted made
by multcomp::glht.

Thank you


-- 
Witold Eryk Wolski

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] convert microns to nm in a messy dataset [solved]

[Ivan Calandra]

Dear Peter,

Thank you for your answer, the function na.locf() is exactly what I needed!
I had started processing my dataset so the first lines (used as headers) were
not included in the sample I have sent. But there is also a "unit" line before
the first value.

And yes, of course, divide by 1000.

Best,
Ivan

--
Dr. Ivan Calandra
TraCEr, laboratory for Traceology and Controlled Experiments
MONREPOS Archaeological Research Centre and
Museum for Human Behavioural Evolution
Schloss Monrepos
56567 Neuwied, Germany
+49 (0) 2631 9772-243
https://www.researchgate.net/profile/Ivan_Calandra

On May 10, 2019 at 3:29 PM peter dalgaard  wrote:
> From nm to micron, _divide_ by 1000 (as you likely know)
>
> What are the units of the first value? Looks like micron in your example, but
> is there a rule?
>
> Basically, it is a "last observation carried forward" type problem, so
> something like this:
>
>
> my.data <- structure(list(V1 = c("2019/05/10", "#", "#", "#", "2019/05/10",
> "2019/05/10", "2019/05/10", "#", "#", "#", "2019/05/10", "#", "#", "#",
> "2019/05/10", "#", "#", "#", "2019/05/10", "2019/05/10"), V19 =
> c("0.2012800083", "45", "Sq", "µm", "0.3634383236", "0.4360454777",
> "0.3767733568", "45", "Sq", "nm", "102.013048", "45", "Sq", "µm",
> "0.1413840498", "45", "Sq", "nm", "65.4459715", "46.45802917")), row.names =
> c(NA, 20L), class = "data.frame")
>
> y <- my.data$V19
> u <- ifelse(y=="nm" | y=="µm", y, NA)
> num <- my.data$V1 != "#"
> uu <- zoo::na.locf(u, na.rm=FALSE)
> data.frame(val = as.numeric(y[num]), units = uu[num])
>
> giving
> val units
> 1 0.2012800 
> 2 0.3634383 µm
> 3 0.4360455 µm
> 4 0.3767734 µm
> 5 102.0130480 nm
> 6 0.1413840 µm
> 7 65.4459715 nm
> 8 46.4580292 nm
>
> and you can surely take it from there.
>
> -pd
>
>
> > On 10 May 2019, at 13:54 , Ivan Calandra  wrote:
> >
> > Dear useRs,
> >
> > Below is a sample of my dataset (I have more rows and columns).
> >
> > As you can see in the 2nd column, there are values, the name of the
> > parameter
> > ('Sq' in that case), some integer ('45' in that case) and the unit ('µm' or
> > 'nm').
> > I know how to extract the rows of interest (those with values), but they are
> > expressed in different units. All values following a line with the unit are
> > expressed in that unit, but the number of lines is not constant (sometimes
> > each
> > value is expressed in a different unit so there will be a new unit line, but
> > there are sometimes several values in a row expressed in the same unit so
> > without unit lines in between). I hope this is clear (it should be with the
> > example provided).
> > This messy dataset comes from an external software so I don't have any means
> > to
> > format the ways the data are collated. I have to find a way to deal with it
> > in
> > R.
> >
> > What I would like to do is convert the values in nm to µm; I just need to
> > multiply by 1000.
> >
> > What I don't know is how to identify the values that are expressed in nm
> > (all
> > values that follow a line with 'nm' until there is a line with 'µm').
> >
> > I don't even know how I should search online because I don't know how this
> > kind
> > of operation is called.
> > Any help is appreciated.
> >
> > Thank you in advance.
> > Ivan
> >
> >
> > my.data <- structure(list(V1 = c("2019/05/10", "#", "#", "#", "2019/05/10",
> > "2019/05/10", "2019/05/10", "#", "#", "#", "2019/05/10", "#", "#", "#",
> > "2019/05/10", "#", "#", "#", "2019/05/10", "2019/05/10"), V19 =
> > c("0.2012800083", "45", "Sq", "µm", "0.3634383236", "0.4360454777",
> > "0.3767733568", "45", "Sq", "nm", "102.013048", "45", "Sq", "µm",
> > "0.1413840498", "45", "Sq", "nm", "65.4459715", "46.45802917")), row.names =
> > c(NA, 20L), class = "data.frame")
> >
> > --
> > Dr. Ivan Calandra
> > TraCEr, laboratory for Traceology and Controlled Experiments
> > MONREPOS Archaeological Research Centre and
> > Museum for Human Behavioural Evolution
> > Schloss Monrepos
> > 56567 Neuwied, Germany
> > +49 (0) 2631 9772-243
> > https://www.researchgate.net/profile/Ivan_Calandra
> >
> > __
> > R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> > 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.
>
> --
> Peter Dalgaard, Professor,
> Center for Statistics, Copenhagen Business School
> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
> Phone: (+45)38153501
> Office: A 4.23
> Email: pd@cbs.dk Priv: pda...@gmail.com
>
>
>
>
>
>
>
>
>
[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] convert microns to nm in a messy dataset

[Ivan Calandra]

Dear John,

Thank you for your answer.
However, it does not make sense to me, as it works only line by line of the
data.frame, and I need something for "last observation carried forward" as Peter
mentioned. The script does not work as is either, probably due to typos with
semi-colons and "if... else" statements, so I cannot really test it.

Best,
Ivan

--
Dr. Ivan Calandra
TraCEr, laboratory for Traceology and Controlled Experiments
MONREPOS Archaeological Research Centre and
Museum for Human Behavioural Evolution
Schloss Monrepos
56567 Neuwied, Germany
+49 (0) 2631 9772-243
https://www.researchgate.net/profile/Ivan_Calandra

On May 10, 2019 at 3:47 PM John McKown  wrote:

>  This is my approach. It is based entirely on what you said (multiply by 1000
> to convert from nm to  µm, but I think it is divide by). It assumes that the
> starting value is in  µm. If the starting value is in nm, change the
> "factor=1" to "factor=1000". That  µm is micro-meters (10^-6) and nm is
> nano-meters (10^-9), so divide by would be correct.
> 
>   factor=1;
>   for (i in 1:length(my.data$V19)) {
>   print("Start");print(factor);print(my.data[i,]);
>   if (my.data$V19[i] == "nm") { factor=1000; my.data$V19[i]="µm";print("nm");}
> else if (my.data$V19[i] == "µm") {factor=1;};
>   if (suppressWarnings(! is.na  (as.numeric(my.data$V19[i] {
> my.data$V19[i] = as.character(as.numeric(my.data$V19[i]) * factor);
> print("changed"); }
>   print(factor);print(my.data[i,]);print("End");
> 
> 
>
>  On Fri, May 10, 2019 at 6:54 AM Ivan Calandra < calan...@rgzm.de
>  > wrote:
>> > Dear useRs,
> >
> >Below is a sample of my dataset (I have more rows and columns).
> >
> >As you can see in the 2nd column, there are values, the name of the
> > parameter
> >('Sq' in that case), some integer ('45' in that case) and the unit ('µm'
> > or
> >'nm').
> >I know how to extract the rows of interest (those with values), but they
> > are
> >expressed in different units. All values following a line with the unit
> > are
> >expressed in that unit, but the number of lines is not constant
> > (sometimes each
> >value is expressed in a different unit so there will be a new unit line,
> > but
> >there are sometimes several values in a row expressed in the same unit so
> >without unit lines in between). I hope this is clear (it should be with
> > the
> >example provided).
> >This messy dataset comes from an external software so I don't have any
> > means to
> >format the ways the data are collated. I have to find a way to deal with
> > it in
> >R.
> >
> >What I would like to do is convert the values in nm to µm; I just need to
> >multiply by 1000.
> >
> >What I don't know is how to identify the values that are expressed in nm
> > (all
> >values that follow a line with 'nm' until there is a line with 'µm').
> >
> >I don't even know how I should search online because I don't know how
> > this kind
> >of operation is called.
> >Any help is appreciated.
> >
> >Thank you in advance.
> >Ivan
> >
> >
> >my.data <- structure(list(V1 = c("2019/05/10", "#", "#", "#",
> > "2019/05/10",
> >"2019/05/10", "2019/05/10", "#", "#", "#", "2019/05/10", "#", "#", "#",
> >"2019/05/10", "#", "#", "#", "2019/05/10", "2019/05/10"), V19 =
> >c("0.2012800083", "45", "Sq", "µm", "0.3634383236", "0.4360454777",
> >"0.3767733568", "45", "Sq", "nm", "102.013048", "45", "Sq", "µm",
> >"0.1413840498", "45", "Sq", "nm", "65.4459715", "46.45802917")),
> > row.names =
> >c(NA, 20L), class = "data.frame")
> >
> >--
> >Dr. Ivan Calandra
> >TraCEr, laboratory for Traceology and Controlled Experiments
> >MONREPOS Archaeological Research Centre and
> >Museum for Human Behavioural Evolution
> >Schloss Monrepos
> >56567 Neuwied, Germany
> >+49 (0) 2631 9772-243
> >https://www.researchgate.net/profile/Ivan_Calandra
> > 
> >
> >__
> >R-help@r-project.org  mailing list -- To
> > UNSUBSCRIBE and more, see
> >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.
> >  >
> 
>  --
>  This is clearly another case of too many mad scientists, and not enough
> hunchbacks.
> 
> 
>  Maranatha! <><
>  John McKown
>

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] multiple graphs on one plot

[Jim Lemon]

Hi Andrew,
First, a little mind reading. My crystal ball says that "cw" can be
interpreted as "carapace width". It didn't tell me the parameters of
the distribution, so:

set.seed(1234)
mf<-list(rnorm(400,145,15),rnorm(400,160,15))
library(plotrix)
multhist(mf, xlab="CW", ylab="Frequency", ylim=c(0,100),main="All Measured
Crabs", col=c("dark gray", "light gray"),
 breaks=seq(90,210, by=10),beside=TRUE,space=c(0,0.5))
legend("topright", c("Females", "Males"), fill=c("dark gray", "light gray"))
lines(seq(0,32,length.out=121),rescale(dnorm(90:210,145,15),c(0,100)))

This produces what I think you are after. Note that it may be
misleading as the distribution of carapace width in real mud crabs
doesn't look normal to me.

Jim

On Mon, May 13, 2019 at 3:00 PM Andrew Halford  wrote:
>
> Hi Listers
>
> I've been trying to make a single graphic that has frequency histograms for
> male and female mud crabs displayed side by side (such as when using the
> beside=TRUE command for barplots). I then want to display a normal
> distribution on top of the male and female histograms.
>
> I have been using the multhist command in Plotrix to generate the
> histograms without too much problem, but I cannot get the normal
> distributions to plot up on the same graph.
>
> Histograms plot
>
> mf <-
> list(lf_crabs$cw[lf_crabs$sex=='female'],lf_crabs$cw[lf_crabs$sex=='male'])
> multhist(mf, xlab="CW", ylab="Frequency", ylim=c(0,100),main="All Measured
> Crabs", col=c("dark gray", "light gray"),
>  breaks=seq(90,210, by=10),beside=TRUE,space=c(0,0.5))
> legend("topright", c("Females", "Males"), fill=c("dark gray", "light gray"))
>
> Then I try to add a normal distribution curve just to the female data but I
> cant get the output to plot
>
> points(seq(min(mf[[1]]), max(mf[[1]]), length.out=300),
>dnorm(seq(min(mf[[1]]), max(mf[[1]]), length.out=300),
>  mean(mf[[1]]), sd(mf[[1]])),type="l", col="dark gray")
>
> Even trying to add an abline to the plot doesn't work.
>
> What am I missing?
>
> cheers
>
> Andy
>
> --
> Andrew Halford Ph.D
> Senior Coastal Fisheries Scientist
> Pacific Community | Communauté du Pacifique CPS – B.P. D5 | 98848 Noumea,
> New Caledonia | Nouméa, Nouvelle-Calédonie
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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 -- To UNSUBSCRIBE and more, see
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.