[R] Integration of functions with a vector argument
Hello all,
I have a question concerning integration of a function of a multivariate
argument with respect to one or more variables in r. Let us say we have a
function
F <- function(x){ body of the function}
Where x is, in general, a d by 1 vector with d>1. Now I want to integrate out
some of the coordinates of x, e.g. x[1] or x[2] or both of them etc. I'm well
aware of how to integrate out e.g. y if a function is defined as f <- function
(x,y) {body of the function} where y is a scalar.
However, it seems to be quite difficult to do the same if the function is
defined with a vector argument x. At the very least, I haven't seen any good
examples of this being done.
Any suggestions?
Yours sincerely,
Michael
Michael Levine
Associate Professor, Statistics
Department of Statistics
Purdue University
250 North University Street
West Lafayette, IN 47907 USA
email: mlev...@purdue.edu
Phone: +1-765-496-7571
Fax: +1-765-494-0558
URL: www.stat.purdue.edu/~mlevins
[[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] R-help Digest, Vol 255, Issue 17
You might be interested in the `Rdatasets` package, https://vincentarelbundock.github.io/Rdatasets/ which lists over 2200 datasets from various packages. What is the context of the `lottery` dataset. I seem to recall smth to do with the NJ Lottery -Michael 1. Availability of Sdatasets (Avro Alo) -- Message: 1 Date: Sun, 19 May 2024 08:58:20 + From: Avro Alo To: "r-help@r-project.org" Subject: [R] Availability of Sdatasets Message-ID: <8I3Bj0m1IzC35J4nEoROCf1yZD66oeLHFLtxsXKSty3vplcl5gKp-_XmdSvEbG0UYtxv8g0Jw0ihsR5x0MS0QdF7DOmooZ2C9BJVqUUlNSQ=@protonmail.com> Content-Type: text/plain; charset="utf-8" >From the mention in R-intro I went to look at The new S language book. In chapter 1 it has a lottery dataset. So naturally I thought it is pre-supplied with R. But I didn't fount, made a google search and found the package that has the dataset, https://docs.tibco.com/pub/enterprise-runtime-for-R/6.1.1/doc/html/Language_Reference/Sdatasets/00Index.html This package is very interesting on it's own. But how can I get it? Also, shouldn't regular R installation have this too? Thanks! (first time posting here) -- Subject: Digest Footer ___ 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. -- End of R-help Digest, Vol 255, Issue 17 __ 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] Strange variable names in factor regression
Comment in in-line below
On 09/05/2024 13:09, Naresh Gurbuxani wrote:
On converting character variables to ordered factors, regression result
has strange names. Is it possible to obtain same variable names with
and without intercept?
Thanks,
Naresh
mydf <- data.frame(date = seq.Date(as.Date("2024-01-01"),
as.Date("2024-03-31"), by = 1))
mydf[, "wday"] <- weekdays(mydf$date, abbreviate = TRUE)
mydf.work <- subset(mydf, !(wday %in% c("Sat", "Sun")))
mydf.weekend <- subset(mydf, wday %in% c("Sat", "Sun"))
mydf.work[, "volume"] <- round(rnorm(nrow(mydf.work), mean = 20, sd =
5))
mydf.weekend[, "volume"] <- round(rnorm(nrow(mydf.weekend), mean = 10,
sd = 5))
mydf <- rbind(mydf.work, mydf.weekend)
reg <- lm(volume ~ wday, data = mydf)
## Variable names as expected
coef(reg)
(Intercept) wdayMon wdaySat wdaySun wdayThu wdayTue
21.3846154 1.3076923 -12.000 -12.9230769 -1.9230769 -0.6923077
wdayWed
-1.6153846
reg <- lm(volume ~ wday - 1, data = mydf)
# Variable names as expected
coef(reg)
wdayFri wdayMon wdaySat wdaySun wdayThu wdayTue wdayWed
21.384615 22.692308 9.384615 8.461538 19.461538 20.692308 19.769231
# Ordered factors for weekday sequence
mydf$wday <- factor(mydf$wday, levels = c("Mon", "Tue", "Wed", "Thu",
"Fri", "Sat", "Sun"), ordered = TRUE)
reg <- lm(volume ~ wday - 1, data = mydf)
# Variable names as expected
coef(reg)
wdayMon wdayTue wdayWed wdayThu wdayFri wdaySat wdaySun
22.692308 20.692308 19.769231 19.461538 21.384615 9.384615 8.461538
reg <- lm(volume ~ wday, data = mydf)
# Strange variable names
coef(reg)
(Intercept) wday.L wday.Q wday.C wday^4 wday^5
17.406593 -12.036715 -4.968654 -1.852819 3.291477 4.263642
wday^6
2.591317
Yes, that is what ordered is supposed to do, fit polynomial contrasts.
When you remove the intercept that breaks it so you get a different fit,
in fact the same as you got when it was not ordered.
Michael
__
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.
--
Michael
__
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] Double buffering plots on Windows
Hi Paul Is there a concrete working example somewhere that shows how to use these to do an animation on Windows (R Gui &/or RStudio) using base R plot() and friends? I have several old examples somewhere that used to work (R < ~ 3), but now no longer work as before. Date: Mon, 25 Mar 2024 10:43:29 +1300 From: Paul Murrell mailto:p...@stat.auckland.ac.nz>> To: "Bickis, Mikelis" mailto:bic...@math.usask.ca>>, "r-help@r-project.org<mailto:r-help@r-project.org>" mailto:r-help@r-project.org>> Subject: Re: [R] Double buffering plots on Windows Message-ID: mailto:b74c68da-a0b2-47dd-b54f-6b318488c...@stat.auckland.ac.nz>> Content-Type: text/plain; charset="utf-8"; Format="flowed" Hi Take a look at dev.hold() and dev.flush() Paul --- Michael Friendly Email: friendly AT yorku DOT ca Professor, Psychology Dept. York University Voice: 416 736-2100 x66249 4700 Keele StreetWeb: http://www.datavis.ca<http://www.datavis.ca/> | @datavisFriendly Toronto, ONT M3J 1P3 CANADA [[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: [ESS] FW: [GNU ELPA] ESS version 24.1.1
On Thu, Feb 22, 2024 at 04:27:18PM +, Sparapani, Rodney via ESS-help wrote: Hello, Rodney! > I see this function in etc/ESSR/R/.load.R on line 12. > So could you please double-check? Thanks I don't see this file in the installed directory, /usr/share/emacs/etc/ess/ESSR/R, but it is in the source directory, /usr/src/ESS-24.01.1/etc/ESSR/R/.load.R as you said, line 12. If I remember rightly, the installation script copied the files to the destination site above. These two .load.R files on my system appear to be very different! I installed ESS from the source with "make" and "make install". Should I overwrite the one with the other? Cheers, Mike -- Dr. Michael L. Dowling Gaußstr. 27 38106 Braunschweig Germany __ ESS-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/ess-help
Re: [ESS] FW: [GNU ELPA] ESS version 24.1.1
On Thu, Feb 15, 2024 at 04:19:43PM +, Sparapani, Rodney via ESS-help wrote: > Hi Gang: > > We just tidied up a recent insidious bug. No other changes. > > Version 24.1.1 of package ESS has just been released in GNU ELPA. > You can now find it in M-x list-packages RET. I have just installed ESS-24.01.1, and it appears to be working. At least, when I run an older program, I get the desired results. But I do get an error message in starting R. It appears only briefly, and disappears immediately when I try to cut it out with cut-and-paste using the mouse. The error message is: Error: could not find function \".ess.ESSR.load\" Emacs appears to be looking for this strangely names file in the directory, /usr/share/emacs/etc/ess/ESSR/R where it clearly is not to be found. Nor can I find any reference to such a function in the source directory, /usr/src/ESS-24.01.1/etc/ESSR/R Can anybody suggest the cause? Cheers, M. Dowling -- Dr. Michael L. Dowling Gaußstr. 27 38106 Braunschweig Germany __ ESS-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/ess-help
[R] DescTools::Quantile
Greetings,
I am having a problem with DescTools::Quantile
(a function computing quantiles from weighted samples):
# these sum to one
probWeights = c(
0.0043, 0.0062, 0.0087, 0.0119, 0.0157, 0.0204, 0.0257, 0.0315, 0.0378,
0.0441, 0.0501, 0.0556, 0.06, 0.0632, 0.0648, 0.0648, 0.0632, 0.06,
0.0556, 0.0501, 0.0441, 0.0378, 0.0315, 0.0257, 0.0204, 0.0157, 0.0119,
0.0087, 0.0062, 0.0043
)
x = seq(-100,100,length.out=length(probWeights))
qtls <- DescTools::Quantile(x, weights=probWeights, probs=c(0.1,0.9))
cat("\nQuantiles:\n")
print(qtls)
Both quantiles are equal to 100!
Is this function working or am I not using it correctly?
Michael Meyer
__
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] Quantiles of sums of independent discrete random variables
Greetings,
Minimal reproducible example as requested by the technical expert Jeff
Newmiller:
library(bayesmeta)
# density of $(1/10)*\sum_{j=1}{10}N(j,0.01$
# (convex sum of normal distributions)
#
f <- Vectorize(function(s) sum(vapply(1:10,
FUN = function(j) dnorm(s,mean=j,sd=0.01)/10, FUN.VALUE=0
)))
g <- function(s) dnorm(s,mean=0,sd=0.01)
cat("\n\n")
for(i in 1:5){
cat("Doing convolution ",i,"\n")
g <- convolve(g,f)$density
}
cat("\nConvolutions finished, plotting density.")
s <- seq(0,100,length.out=1024)
matplot(s,g(s),type="l")
Michael Meyer
[[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] Use of geometric mean .. in good data analysis
By the Strong Law of Large Numbers applied to log(X) the geometric mean of X_1,...,X_n > 0 and IID like X converges toexp(E[log(X)]] which, by Jensen's inequality, is always <= E[X] and is strictly less than E[X] except in trivial extreme cases. In short: by using the geometric mean all asymptotic results no longer apply. Michael Meyer [[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] Quantiles of sums of independent discrete random variables
Greetings, I have the following Problem: Given k (=10) discrete independent random variables X_i with n_i (= 5 to 20) values each,compute quantiles of the distribution of the sum X = X_1+...+X_k. Here X has n=n_1 x n_2 ... n_k distinct values which is too large to list them all together with their probabilities. I tried several approaches: (A) Convolution: each X_j is approximated with Y_j=X_j+Z, where Z is an N(0,sigma) variable with small sigma. Then Y_j is a probability mixture of the normal variables N(x_j,sigma), where the x_j runs over all values of X_j, and has a highly oscillatory density. The density of Y=\sum Y_j is the convolution of the densities of the Y_j. I need this density at a sizeable number of points and this turns out to be too slow. The issue seems to be the convergence of the convolution integrals slowed down by the oscillatory nature of the densities of the Y_j. When the densities are behaved better (e.g. normal RVs), the computation of such a convolution is quite fast. (B) Characteristic function: X will be approximated with Y=X+Z, where Z is normal N(0,sigma) with small sigma. Y has a density (which it is impossible to compute directly) but the characteristic function (_continuous_ Fourier transform) cf_Y of Y can easily be computed analytically (without knowing the density of Y) Now let s be a numeric vector. I want to get the density f_Y(s) of Y evaluated along s. The proper way of doing this would be to apply the inverse continuous Fourier transform to the function cf_Y at each point in s. This is far too slow. That's why I tried to apply the inverse _discrete_ Fourier transform to the vector of values cf_Y(s) and that does not yield anything reasonable. This baffles me since I was under the impression that the discrete Fourier transform is an approximation to the continuous Fourier transform and so should yield the values of the latter, if the value grid s is fine enough. Should this work? Note that this inversion would definitely work if I could compute the _discrete_ Fourier transform of the density f_Y along s, but regrettably this is not possible since (a) the density of f_Y is far too complicated, and (b) the discrete Fourier transform of a sum Y = Y_1+Y_2+...+Y_k of independent random variables Y_j is not the product of the discrete Fourier transforms of the Y_j. Any ideas how I could approach this problem with the tools of R? Thanks in advance for all replies, Michael Meyer [[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] Function with large nested list
Dear Emily
Comment in-line
On 18/12/2023 09:56, Emily Bakker wrote:
Hello list,
I want to make a large rulebased algorithm, to provide decision support for
drug prescriptions. I have defined the algorithm in a function, with a for loop
and many if statements. The structure should be as follows:
1. Iterate over a list of drug names. For each drug:
2. Get some drug related data (external dataset). Row of a dataframe.
3. Check if adaptions should be made to standard dosage and safety information
in case of contraindications. If patient has an indication, update current
dosage and safety information with the value from the dataframe row.
4. Save dosage and safety information in some lists and continue to the next
drug.
5. When the iteration over all drugs is done, return the lists.
ISSUE:
So it is a very large function with many nested if statements. I have checked the code structure multiple times,
but i run into some issues. When i try to run the function definiton, the command never "completes" in
de console. Instead of ">", the console shows "+". No errors are raised.
When my console returns a + is usually means I have left off the final
parenthesis or given it an incomplete line.
Michael
As I said, i have checked the structure multiple times, but cant find an error.
I have tried rebuilding it and testing each time i add a part. Each part
functions isolated, but not together in the same function. I can't find any
infinite loops either.
I suspect the function may be too large, and i have to define functions for
each part separately. That isn't an issue necessarily, but i would still like
to know why my code won't run. And whether there are any downsides or
considerations for using many small functions.
Below is my code. I have left part of it out. There are six more parts like the
diabetes part that are similar.
I also use a lot of data/variabeles not included here, to try and keep things
compact. But I can provide additional information if helpful.
Thanks it advance for thinking along!!
Kind regards,
Emily
The code:
decision_algorithm <- function(AB_list, dataset_ab = data.frame(), diagnose = 'cystitis',
diabetes_status = "nee", katheter_status = "nee",
lang_QT_status = "nee", obesitas_status = "nee",
zwangerschap_status = "nee",
medicatie_actief = data.frame(dict[["med_AB"]]), geslacht
= "man", gfr=90){
# vars
list_AB_status <- setNames(as.list(rep("green", length(AB_list))),
names(AB_list)) #make a dict of all AB's and assign status green as deafault for status
list_AB_remarks <- setNames(as.list(rep("Geen opmerkingen", length(AB_list))),
names(AB_list)) #make a dict of all AB's and assign "Geen" as default for remarks #Try empty
list
list_AB_dosering <- setNames(as.list(rep("Geen informatie", length(AB_list))),
names(AB_list)) # make named list of all AB's and assign "Geen informatie", will be replaced
with actual information in algorithm
list_AB_duur <- setNames(as.list(rep("Geen informatie", length(AB_list))),
names(AB_list)) # make named list of all AB's and assign "Geen informatie", will be replaced
with actual information in algorithm
# CULTURES #
for (i in names(AB_list)) {
ab_data <- dataset_ab[dataset_ab$middel == i,] #get info for this AB from dataset_ab
# Extract and split the diagnoses, dosering, and duur info for the current antibiotic
ab_diagnoses <- str_split(ab_data$diagnoses, pattern = " \\| ")[[1]]
ab_diagnose_dosering <- str_split(ab_data$`diagnose dosering`, pattern = " \\|
")[[1]]
ab_diagnose_duur <- str_split(ab_data$`diagnose duur`, pattern = " \\|
")[[1]]
# Find the index of the current diagnose in the ab_diagnoses list
diagnose_index <- match(diagnose, ab_diagnoses)
# Determine dosering and duur based on the diagnose_index
if (!is.na(diagnose_index)) {
dosering <- ifelse(ab_diagnose_dosering[diagnose_index] == "standaard",
ab_data$dosering, ab_diagnose_dosering[diagnose_index])
duur <- ifelse(ab_diagnose_duur[diagnose_index] == "standaard",
ab_data$duur, ab_diagnose_duur[diagnose_index])
} else {
# Use general dosering and duur as fallback if diagnose is not found
dosering <- ab_data$dosering
duur <- ab_data$duur
}
list_AB_dosering[[i]] <- dosering
list_AB_duur[[i]] <- duur
if ((!is.null(AB_list[[i]]) && AB_list[[i]] == "I")) {
list_AB_status[[i]] <- "yellow"
list_AB_remarks[[i]] <- "Kweek verminderd gevoelig"
} else if ((!is.n
Re: [R] [External] Re: Error in setwd(dir) when initializing R
Dear Ana When you installed R it gave you an option to put a shortcut on the desktop. If you did that then you can start the R GUI by clicking on it. You may need to customise it for your preferences. Specifically it starts R with an option cd-user-docs which I delete and I then alter the starting option to start in the directory I want. You may have different preferences of course. Try Rgui as Ivan suggests and report back what happened. You may have to ask your local IT support as if this is something the Complutense is setting up it is surprising that all your collegues are free of problems. Have you asked around? Michael On 21/11/2023 12:41, Ana de las Heras Molina wrote: Hello, I am sorry for my ignorance, but what is Rgui.exe? El mar, 21 nov 2023 a las 12:06, Ivan Krylov () escribió: В Tue, 21 Nov 2023 10:51:59 +0100 Ana de las Heras Molina пишет: I uninstalled onedrive, I eliminated all the folders and then reinstalled R and RStudio... but it is RStudio the one creating a folder called C:\Users\Ana\OneDrive - Universidad Complutense de Madrid (UCM)\Documentos. Does Rgui.exe use a different home directory from RStudio? Perhaps you need not to reinstall RStudio (which deletes and recreates the program files which aren't damaged) but to wipe the profile (the settings somewhere under %APPDATA% or %LOCALAPPDATA%), but it pays to be extra careful about these directories. People at <https://community.rstudio.com/> may know more about that. I'm afraid I don't know how to disable OneDrive on a Windows 10 installation. If your university has a Microsoft support contract, it could be worth asking the Microsoft representative to fix OneDrive so that temporary files would work on it and telling them the steps to reproduce the problem. Browse[1]> dir.exists(dir) [1] TRUE Browse[1]> setwd(dir) Error durante el wrapup: no es posible cambiar el directorio de trabajo What is the value of `dir` at this point? (What does print(dir) say?) Can you open this directory in Windows Explorer? If possible, could you please set your mailer to compose messages in plain text instead of HTML? The plain text versions of your messages that your mailer generates from the HTML messages are somewhat mangled: https://stat.ethz.ch/pipermail/r-help/2023-November/478593.html -- Best regards, Ivan -- Michael __ 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] car::deltaMethod() fails when a particular combination of categorical variables is not present
Thank you very much, John. This has allowed us to move forward on several fronts and better understand our data. - Michael Cohn On Tue, Sep 26, 2023 at 8:39 AM John Fox wrote: > Dear Michael, > > My previous response was inaccurate: First, linearHypothesis() *is* able > to accommodate aliased coefficients by setting the argument singular.ok > = TRUE: > > > linearHypothesis(minimal_model, "bt2 + csent + bt2:csent = 0", > + singular.ok=TRUE) > > Linear hypothesis test: > bt2 + csent + bt2:csent = 0 > > Model 1: restricted model > Model 2: a ~ b * c > >Res.DfRSS Df Sum of Sq F Pr(>F) > 1 16 9392.1 > 2 15 9266.4 1125.67 0.2034 0.6584 > > Moreover, when there is an empty cell, this F-test is (for a reason that > I haven't worked out, but is almost surely due to how the rank-deficient > model is parametrized) *not* equivalent to the t-test for the > corresponding coefficient in the raveled version of the two factors: > > > df$bc <- factor(with(df, paste(b, c, sep=":"))) > > m <- lm(a ~ bc, data=df) > > summary(m) > > Call: > lm(formula = a ~ bc, data = df) > > Residuals: > Min 1Q Median 3Q Max > -57.455 -11.750 0.439 14.011 37.545 > > Coefficients: > Estimate Std. Error t value Pr(>|t|) > (Intercept)20.50 17.57 1.166 0.2617 > bct1:unsent37.50 24.85 1.509 0.1521 > bct2:other 32.00 24.85 1.287 0.2174 > bct2:sent 17.17 22.69 0.757 0.4610 <<< cf. F = 0.2034, p > = 0.6584 > bct2:unsent38.95 19.11 2.039 0.0595 > > Residual standard error: 24.85 on 15 degrees of freedom > Multiple R-squared: 0.2613,Adjusted R-squared: 0.06437 > F-statistic: 1.327 on 4 and 15 DF, p-value: 0.3052 > > In the full-rank case, however, what I said is correct -- that is, the > F-test for the 1 df hypothesis on the three coefficients is equivalent > to the t-test for the corresponding coefficient when the two factors are > raveled: > > > linearHypothesis(minimal_model_fixed, "bt2 + csent + bt2:csent = 0") > > Linear hypothesis test: > bt2 + csent + bt2:csent = 0 > > Model 1: restricted model > Model 2: a ~ b * c > >Res.DfRSS Df Sum of Sq F Pr(>F) > 1 15 9714.5 > 2 14 9194.4 1520.08 0.7919 0.3886 > > > df_fixed$bc <- factor(with(df_fixed, paste(b, c, sep=":"))) > > m <- lm(a ~ bc, data=df_fixed) > > summary(m) > > Call: > lm(formula = a ~ bc, data = df_fixed) > > Residuals: > Min 1Q Median 3Q Max > -57.455 -11.750 0.167 14.011 37.545 > > Coefficients: > Estimate Std. Error t value Pr(>|t|) > (Intercept) 64.000 25.627 2.497 0.0256 > bct1:sent-43.500 31.387 -1.386 0.1874 > bct1:unsent -12.000 36.242 -0.331 0.7455 > bct2:other -11.500 31.387 -0.366 0.7195 > bct2:sent-26.333 29.591 -0.890 0.3886 << cf. > bct2:unsent -4.545 26.767 -0.170 0.8676 > > Residual standard error: 25.63 on 14 degrees of freedom > Multiple R-squared: 0.2671,Adjusted R-squared: 0.005328 > F-statistic: 1.02 on 5 and 14 DF, p-value: 0.4425 > > So, to summarize: > > (1) You can use linearHypothesis() with singular.ok=TRUE to test the > hypothesis that you specified, though I suspect that this hypothesis > probably isn't testing what you think in the rank-deficient case. I > suspect that the hypothesis that you want to test is obtained by > raveling the two factors. > > (2) There is no reason to use deltaMethod() for a linear hypothesis, but > there is also no intrinsic reason that deltaMethod() shouldn't be able > to handle a rank-deficient model. We'll probably fix that. > > My apologies for the confusion, > John > > -- > John Fox, Professor Emeritus > McMaster University > Hamilton, Ontario, Canada > web: https://www.john-fox.ca/ > > On 2023-09-26 9:49 a.m., John Fox wrote: > > Caution: External email. > > > > > > Dear Michael, > > > > You're testing a linear hypothesis, so there's no need to use the delta > > method, but the linearHypothesis() function in the car package also > > fails in your case: > > > > > linearHypothesis(minimal_model, "bt2 + csent + bt2:csent = 0") > > Error in linearHypothesis.lm(minimal_model, "bt2 + csent + bt2:csent = > > 0") : > > there are aliased coefficients in the model. > > > > One work-around is to ravel the two factors into a single factor with 5 > > levels: > > > > > df$bc <- fact
Re: [R] Question about R software and output
Dear Charity Since your organisation is a member of King's Health Partners you might like to ask colleagues in KCL for local support. Michael On 02/10/2023 08:48, Ferguson Charity (CEMINFERGUSON) wrote: To whom it may concern, My understanding is that the R software is downloaded from a CRAN network and data is imported into it using Microsoft Excel for example. Could I please just double check whether any data or results from the output is held on external servers or is it just held on local files on the computer? Many thanks, Charity * The information contained in this message and or attachments is intended only for the person or entity to which it is addressed and may contain confidential and/or privileged material. Unless otherwise specified, the opinions expressed herein do not necessarily represent those of Guy's and St Thomas' NHS Foundation Trust or any of its subsidiaries. The information contained in this e-mail may be subject to public disclosure under the Freedom of Information Act 2000. Unless the information is legally exempt from disclosure, the confidentiality of this e-mail and any replies cannot be guaranteed. Any review, retransmission,dissemination or other use of, or taking of any action in reliance upon, this information by persons or entities other than the intended recipient is prohibited. If you received this in error, please contact the sender and delete the material from any system and destroy any copies. We make every effort to keep our network free from viruses. However, it is your responsibility to ensure that this e-mail and any attachments are free of viruses as we can take no responsibility for any computer virus which might be transferred by way of this e-mail. * [[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. -- Michael __ 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] car::deltaMethod() fails when a particular combination of categorical variables is not present
I'm running a linear regression with two categorical predictors and their
interaction. One combination of levels does not occur in the data, and as
expected, no parameter is estimated for it. I now want to significance test
a particular combination of levels that does occur in the data (ie, I want
to get a confidence interval for the total prediction at given levels of
each variable).
In the past I've done this using car::deltaMethod() but in this dataset
that does not work, as shown in the example below: The regression model
gives the expected output, but deltaMethod() gives this error:
error in t(gd) %*% vcov. : non-conformable arguments
I believe this is because there is no parameter estimate for when the
predictors have the values 't1' and 'other'. In the df_fixed dataframe,
putting one person into that combination of categories causes deltaMethod()
to work as expected.
I don't know of any theoretical reason that missing one interaction
parameter estimate should prevent getting a confidence interval for a
different combination of predictors. Is there a way to use deltaMethod() or
some other function to do this without changing my data?
Thank you,
- Michael Cohn
Vote Rev (http://voterev.org)
Demonstration:
--
library(car)
# create dataset with outcome and two categorical predictors
outcomes <- c(91,2,60,53,38,78,48,33,97,41,64,84,64,8,66,41,52,18,57,34)
persontype <-
c("t2","t2","t2","t2","t2","t2","t2","t1","t2","t2","t2","t2","t1","t1","t2","t2","t1","t2","t2","t2")
arm_letter <-
c("unsent","unsent","unsent","unsent","sent","unsent","unsent","sent","unsent","unsent","other","unsent","unsent","sent","unsent","other","unsent","sent","sent","unsent")
df <- data.frame(a = outcomes, b=persontype, c=arm_letter)
# note: there are no records with the combination 't1' + 'other'
table(df$b,df$c)
#regression works as expected
minimal_formula <- formula("a ~ b*c")
minimal_model <- lm(minimal_formula, data=df)
summary(minimal_model)
#use deltaMethod() to get a prediction for individuals with the combination
'b2' and 'sent'
# deltaMethod() fails with "error in t(gd) %*% vcov. : non-conformable
arguments."
deltaMethod(minimal_model, "bt2 + csent + `bt2:csent`", rhs=0)
# duplicate the dataset and change one record to be in the previously empty
cell
df_fixed <- df
df_fixed[c(13),"c"] <- 'other'
table(df_fixed$b,df_fixed$c)
#deltaMethod() now works
minimal_model_fixed <- lm(minimal_formula, data=df_fixed)
deltaMethod(minimal_model_fixed, "bt2 + csent + `bt2:csent`", rhs=0)
[[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] How to fix this problem
It looks here as though the E coli column has commas in it so will be treated as character. Michael On 25/09/2023 15:45, avi.e.gr...@gmail.com wrote: David, This may just be the same as your earlier problem. When the type of a column is guessed by looking at the early entries, any non-numeric entry forces the entire column to be character. Suggestion: fix your original EXCEL FILE or edit your CSV to remove the last entries that look just lie commas. -Original Message- From: R-help On Behalf Of Parkhurst, David Sent: Sunday, September 24, 2023 2:06 PM To: r-help@r-project.org Subject: [R] How to fix this problem I have a matrix, KD6, and I�m trying to get a correlation matrix from it. When I enter cor(KD6), I get the message �Error in cor(KD6) : 'x' must be numeric�. Here are some early lines from KD6: Flow E..coliTNSRP TPTSS 1 38.82,4201.65300 0.0270 0.0630 66.80 2 133.02,4201.39400 0.0670 0.1360 6.80 3 86.2 101.73400 0.0700 0.1720 97.30 4 4.85,3900.40400 0.0060 0.0280 8.50 5 0.32,4900.45800 0.0050 0.0430 19.75 6 0.0 1860.51200 0.0040 0.0470 12.00 7 11.19,8351.25500 0.0660 0.1450 12.20 Why are these not numeric? There are some NAs later in the matrix, but I get this same error if I ask for cor(KD6[1:39,]) to leave out the lines with NAs. Are they a problem anyway? [[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. -- Michael __ 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] Odd result
Dear David To get the first 46 rows just do KurtzData[1:43,] However really you want to find out why it happened. It looks as though the .csv file you read has lots of blank lines at the end. I would open it in an editor to check that. Michael On 23/09/2023 23:55, Parkhurst, David wrote: With help from several people, I used file.choose() to get my file name, and read.csv() to read in the file as KurtzData. Then when I print KurtzData, the last several lines look like this: 39 5/31/22 16.0 3411.75525 0.0201 0.0214 7.00 40 6/28/22 2:00 PM 0.0 2150.67950 0.0156 0.0294 NA 41 7/25/22 11:00 AM 11.9 1943.5NA NA 0.0500 7.80 42 8/31/22 0220.5NA NA 0.0700 30.50 43 9/28/22 0.067 10.9NA NA 0.0700 10.20 44 10/26/22 0.086 237NA NA 0.1550 45.00 45 1/12/23 1:00 PM 36.2624196NA NA 0.7500 283.50 46 2/14/23 1:00 PM 20.71 55NA NA 0.0500 2.40 47 NA NA NA NA 48 NA NA NA NA 49 NA NA NA NA Then the NA�s go down to one numbered 973. Where did those extras likely come from, and how do I get rid of them? I assume I need to get rid of all the lines after #46, to do calculations and graphics, no? David [[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. -- Michael __ 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] prop.trend.test
Dear Thomas Are you looking for more than smokers / patients Michael On 07/09/2023 14:23, Thomas Subia via R-help wrote: Colleagues Consider smokers <- c( 83, 90, 129, 70 ) patients <- c( 86, 93, 136, 82 ) prop.trend.test(smokers, patients) Output: Chi-squared Test for Trend inProportions data: smokers out of patients , using scores: 1 2 3 4 X-squared = 8.2249, df = 1, p-value = 0.004132 # trend test for proportions indicates proportions aretrending. How does one identify the direction of trending? # prop.test indicates that the proportions are unequal but doeslittle to indicate trend direction. All the best, Thomas Subia [[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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] To UNSUBSCRIBE
You can find instructions on how to do this on the bottom of e-mailas from the list. On 25/06/2023 15:34, Siddharth Arun wrote: Please unsubscribe me from your mailing list. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] post tweets with media attachment - and with retirement of rtweet
All - I run an automated script that pulls in some market data, makes some charts, and posts them to twitter. A couple of them per day, and then once per month there's a little flurry around the CPI report. For a while rtweet() worked great. Then v2 happened and Oauth2.0. The package broke, was fixed, and then broke (on 1.1 endpoints) for good. It works with 2.0 endpoints, but I can only post text...not media. Most of the R packages for Twitter focus on retrieving tweets and information from Twitter, and not on posting to Twitter. So the options were already pretty limited. And now, even if I were to get rtweet to work for me (with respect to attaching media...and is there any way to auto-refresh the Oauth2.0 token without manually re-authorizing it every two hours???) it appears that rtweet is no longer going to be updated. The immediate issue is figuring out how to get rtweet to work for me to post media along with the text post. But maybe there's a better package. Can anyone suggest a course of action? (I guess worst case is that I continue to generate the media with R, then switch over to Python or something to post the tweet, but that's far more manual than I want - the whole point is for this to run unattended.) Open to ideas! I've been wrestling with this for three days and finally realized that r-help is a resource. All ideas welcome! Mike [[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] Error on using Confint function to calculate 95% CI
I am afraid your post is more or less unreadable since you posted in HTML and this is a plain text list. It might also help if you gave more context like the full results of your model. There is a list dedicated to mixed models https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models which may be able to help you better. Michael On 26/04/2023 02:37, bharat rawlley via R-help wrote: Hello, I am trying to estimate 95% CI using the confint function for a generalized liner model. While I am able to estimate Odds ratio using the coef function but on using the confint function, I get the message " Error in approx(sp$y, sp$x, xout = cutoff) : need at least two non-NA values to interpolate" The following are the Coefficients for which I am trying to calculate 95% CI. The original data has no NA values so I am not sure why this error is showing up Estimate Std. Error z value Pr(>|z|)1.05649 0.30658 3.446 0.000569 ***-0.58348 0.49948 -1.168 0.242744-0.77959 1.47292 -0.529 0.596609-18.50414 6522.63862 -0.003 0.997736-0.82328 0.49102 -1.677 0.093608 .-13.95852 4036.00271 -0.003 0.9972410.53909 0.83073 0.649 0.516376-0.55514 0.46189 -1.202 0.22940817.75434 6522.63863 0.003 0.997828-15.97188 1051.20492 -0.015 0.987878-0.16589 0.43172 -0.384 0.70079231.60387 1459.11421 0.022 0.982719-0.42616 0.43637 -0.977 0.3287661.10179 1.19697 0.920 0.3573210.52345 1.23113 0.425 0.6707050.50791 1.54670 0.328 0.7426210.19881 0.34945 0.569 0.569403-0.07887 0.49769 -0.158 0.8740850.05958 0.36780 0.162 0.871306-0.61339 0.38503 -1.593 0.42 I tried to compute 95% CI using the above coefficients using the cofit function Thank you! [[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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] nlme varFixed
Dear Franz Only attachments of specific types are allowed on the R lists and yours did not come through. It might be necessary to show your code as well. Have you thought of trying the mailing list dedicated to mixed models? https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models Michael On 04/03/2023 09:16, Franz Härtl wrote: Dear R-project team, I have a problem with the function varFixed() of the nlme-package. I used it with the squid-data of Zuur et. al 2009 (chapter 4), to fix increasing residuals (heterogenetiy) (see graph in the email) I get the message ' Variance function structure of class varFixed with no parameters, or uninitialized Could you help me please? Kind regards Franz __ 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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] dyn.load(now = FALSE) not actually lazy?
g
file=/stornext/System/data/apps/pcre2/pcre2-10.36-gcc-11.1.0/lib/libplc4.so
1655: trying
file=/stornext/System/data/tools/openSSL/openSSL-1.1.1n/lib/libplc4.so
1655: trying
file=/stornext/System/data/apps/R/R-4.2.1/lib64/R/lib/libplc4.so
1655: trying file=/usr/local/lib64/libplc4.so
1655: trying
file=/stornext/System/data/tools/openjdk/openjdk-13.0.2/lib/server/libplc4.so
1655: trying
file=/stornext/System/data/apps/hdf5/hdf5-1.8.20/lib/libplc4.so
1655: search cache=/etc/ld.so.cache
1655: trying file=/lib64/libplc4.so
1655:
1655: find library=libnspr4.so [0]; searching
1655: search
path=/stornext/System/data/apps/gcc/gcc-11.1.0/lib64:/stornext/System/data/apps/pcre2/pcre2-10.36-gcc-11.1.0/lib:/stornext/System/data/tools/openSSL/openSSL-1.1.1n/lib:/stornext/System/data/apps/R/R-4.2.1/lib64/R/lib:/usr/local/lib64:/stornext/System/data/tools/openjdk/openjdk-13.0.2/lib/server:/stornext/System/data/apps/hdf5/hdf5-1.8.20/lib
(LD_LIBRARY_PATH)
1655: trying
file=/stornext/System/data/apps/gcc/gcc-11.1.0/lib64/libnspr4.so
1655: trying
file=/stornext/System/data/apps/pcre2/pcre2-10.36-gcc-11.1.0/lib/libnspr4.so
1655: trying
file=/stornext/System/data/tools/openSSL/openSSL-1.1.1n/lib/libnspr4.so
1655: trying
file=/stornext/System/data/apps/R/R-4.2.1/lib64/R/lib/libnspr4.so
1655: trying file=/usr/local/lib64/libnspr4.so
1655: trying
file=/stornext/System/data/tools/openjdk/openjdk-13.0.2/lib/server/libnspr4.so
1655: trying
file=/stornext/System/data/apps/hdf5/hdf5-1.8.20/lib/libnspr4.so
1655: search cache=/etc/ld.so.cache
1655: trying file=/lib64/libnspr4.so
1655:
1655: find library=libselinux.so.1 [0]; searching
1655: search
path=/stornext/System/data/apps/gcc/gcc-11.1.0/lib64:/stornext/System/data/apps/pcre2/pcre2-10.36-gcc-11.1.0/lib:/stornext/System/data/tools/openSSL/openSSL-1.1.1n/lib:/stornext/System/data/apps/R/R-4.2.1/lib64/R/lib:/usr/local/lib64:/stornext/System/data/tools/openjdk/openjdk-13.0.2/lib/server:/stornext/System/data/apps/hdf5/hdf5-1.8.20/lib
(LD_LIBRARY_PATH)
1655: trying
file=/stornext/System/data/apps/gcc/gcc-11.1.0/lib64/libselinux.so.1
1655: trying
file=/stornext/System/data/apps/pcre2/pcre2-10.36-gcc-11.1.0/lib/libselinux.so.1
1655: trying
file=/stornext/System/data/tools/openSSL/openSSL-1.1.1n/lib/libselinux.so.1
1655: trying
file=/stornext/System/data/apps/R/R-4.2.1/lib64/R/lib/libselinux.so.1
1655: trying file=/usr/local/lib64/libselinux.so.1
1655: trying
file=/stornext/System/data/tools/openjdk/openjdk-13.0.2/lib/server/libselinux.so.1
1655: trying
file=/stornext/System/data/apps/hdf5/hdf5-1.8.20/lib/libselinux.so.1
1655: search cache=/etc/ld.so.cache
1655: trying file=/lib64/libselinux.so.1
1655:
1655: find library=libpcre.so.1 [0]; searching
1655: search
path=/stornext/System/data/apps/gcc/gcc-11.1.0/lib64:/stornext/System/data/apps/pcre2/pcre2-10.36-gcc-11.1.0/lib:/stornext/System/data/tools/openSSL/openSSL-1.1.1n/lib:/stornext/System/data/apps/R/R-4.2.1/lib64/R/lib:/usr/local/lib64:/stornext/System/data/tools/openjdk/openjdk-13.0.2/lib/server:/stornext/System/data/apps/hdf5/hdf5-1.8.20/lib
(LD_LIBRARY_PATH)
1655: trying
file=/stornext/System/data/apps/gcc/gcc-11.1.0/lib64/libpcre.so.1
1655: trying
file=/stornext/System/data/apps/pcre2/pcre2-10.36-gcc-11.1.0/lib/libpcre.so.1
1655: trying
file=/stornext/System/data/tools/openSSL/openSSL-1.1.1n/lib/libpcre.so.1
1655: trying
file=/stornext/System/data/apps/R/R-4.2.1/lib64/R/lib/libpcre.so.1
1655: trying file=/usr/local/lib64/libpcre.so.1
1655: trying
file=/stornext/System/data/tools/openjdk/openjdk-13.0.2/lib/server/libpcre.so.1
1655: trying
file=/stornext/System/data/apps/hdf5/hdf5-1.8.20/lib/libpcre.so.1
1655: search cache=/etc/ld.so.cache
1655: trying file=/lib64/libpcre.so.1
1655:
1655:
/stornext/System/data/tools/pgsql/pgsql-15.1/lib/libpqwalreceiver.so:
error: symbol lookup error: undefined symbol: work_mem (fatal)
Error in
dyn.load("/stornext/System/data/tools/pgsql/pgsql-15.1/lib/libpqwalreceiver.so",
:
unable to load shared object
'/stornext/System/data/tools/pgsql/pgsql-15.1/lib/libpqwalreceiver.so':
/stornext/System/data/tools/pgsql/pgsql-15.1/lib/libpqwalreceiver.so:
undefined symbol: work_mem
Execution halted
On Wed, Feb 1, 2023 at 9:05 PM Ivan Krylov wrote:
> В Wed, 1 Feb 2023 14:16:54 +1100
> Michael Milton пишет:
>
> > Is this a bug in the `dyn.load` implementation for R? If not, why is
> > it behaving like this? What should I do about it?
>
> On Unix-like syste
[R] dyn.load(now = FALSE) not actually lazy?
On Linux, if I have a .so file that has a dependency on another .so, and I
`dyn.load(now=FALSE)` the first one, R seems to try to resolve the symbols
immediately, causing the load to fail.
For example, I have `libtorch` installed on my HPC. Note that it links to
various libs such as `libcudart.so` and `libmkl_intel_lp64.so.2` which
aren't currently in my library path:
➜ ~ ldd
/stornext/System/data/nvidia/libtorch-gpu/libtorch-gpu-1.12.1/lib/libtorch_cpu.so
linux-vdso.so.1 => (0x7ffcab58c000)
libgomp.so.1 =>
/stornext/System/data/apps/gcc/gcc-11.2.0/lib64/libgomp.so.1
(0x7f8cb22bf000)
libpthread.so.0 => /lib64/libpthread.so.0 (0x7f8cb20a3000)
libc10.so =>
/stornext/System/data/nvidia/libtorch-gpu/libtorch-gpu-1.12.1/lib/libc10.so
(0x7f8cb1e2d000)
libnuma.so.1 => /lib64/libnuma.so.1 (0x7f8cb1c21000)
librt.so.1 => /lib64/librt.so.1 (0x7f8cb1a19000)
libgcc_s.so.1 =>
/stornext/System/data/apps/gcc/gcc-11.2.0/lib64/libgcc_s.so.1
(0x7f8cb1801000)
libdl.so.2 => /lib64/libdl.so.2 (0x7f8cb15fd000)
libmkl_intel_lp64.so.2 => not found
libmkl_gnu_thread.so.2 => not found
libmkl_core.so.2 => not found
libm.so.6 => /lib64/libm.so.6 (0x7f8cb12fb000)
libcudart.so.11.0 => not found
Then in R, if I try to load that same file:
>
dyn.load("/stornext/System/data/nvidia/libtorch-gpu/libtorch-gpu-1.12.1/lib/libtorch_cpu.so",
now=FALSE)
Error in
dyn.load("/stornext/System/data/nvidia/libtorch-gpu/libtorch-gpu-1.12.1/lib/libtorch_cpu.so",
:
unable to load shared object
'/stornext/System/data/nvidia/libtorch-gpu/libtorch-gpu-1.12.1/lib/libtorch_cpu.so':
libmkl_intel_lp64.so.2: cannot open shared object file: No such file or
directory
Is this a bug in the `dyn.load` implementation for R? If not, why is it
behaving like this? What should I do about it?
For reference, I'm on CentOS 7, with Linux
kernel 3.10.0-1160.81.1.el7.x86_64.
[[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] My forest plot is not fit to windows in R software
Dear Fahimeh
Try preceding your call of forest.meta with
pdf("myplot.pdf")
forest.meta(..)
dev.off()
Why did you set xlim to c(0, 1)
Michael
On 10/01/2023 19:37, Fahimeh Alizadeh wrote:
Dear Prof.
Thank you for your kind response. I have only used:
install.packages("tidyverse")
install.packages("meta")
install.packages("metafor")
library(tidyverse)
library(meta)
library(metafor)
and then,
forest.meta(hidemeta, layout="|RevMan5|", xlab="Proportion", comb.r=T, comb.f=F, xlim = c(0,1), fontsize=10,
digits=3)
Unfortunately I am not familiar to R software, I have not chosen graphics
device. could you please help me to choosegraphics format like pdf?
All the best for you
Sincerely,
Fahimeh Alizadeh
Dr. Fahimeh Alizadeh
Post-Doctorate Fellow and Ph.D
Microbial Biotechnology
Islamic Azad University
On Tuesday, January 10, 2023 at 09:46:40 PM GMT+3:30, Michael Dewey
wrote:
Dear Fahimeh
It would help if you also told us what package you are using as there
are several which perform forest plots. I assume it is meta.
Your plot is cropped on three side as far as I can see. What parameters
did you give to the call of the graphics device? What happens if you use
a different graphics format like pdf? Then we may be able to see where
the issue lies.
Michael
On 10/01/2023 12:30, Fahimeh Alizadeh via R-help wrote:
> Dear Prof. My forest plot is not fit to windows in R software, I have
49 studies but forest plot only show from 9 to 42.
> forest.meta(hidemeta, layout="RevMan5", xlab="Proportion", comb.r=T,
comb.f=F, xlim = c(0,1), fontsize=10, digits=3)
>
>
> Whould you please help me to access full forest plot? Thank you
>
> Sincerely,Fahimeh
>
> Dr. Fahimeh AlizadehPost-Doctorate Fellow and Ph.DMicrobial Biotechnology
> Islamic Azad University
>
>
>
>
>
>
> __
> R-help@r-project.org <mailto:R-help@r-project.org> mailing list -- To
UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
<https://stat.ethz.ch/mailman/listinfo/r-help>
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
<http://www.R-project.org/posting-guide.html>
> and provide commented, minimal, self-contained, reproducible code.
--
Michael
http://www.dewey.myzen.co.uk/home.html
<http://www.dewey.myzen.co.uk/home.html>
<http://www.avg.com/email-signature?utm_medium=email_source=link_campaign=sig-email_content=emailclient>
Virus-free.www.avg.com
<http://www.avg.com/email-signature?utm_medium=email_source=link_campaign=sig-email_content=emailclient>
<#DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2>
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] My forest plot is not fit to windows in R software
Dear Fahimeh It would help if you also told us what package you are using as there are several which perform forest plots. I assume it is meta. Your plot is cropped on three side as far as I can see. What parameters did you give to the call of the graphics device? What happens if you use a different graphics format like pdf? Then we may be able to see where the issue lies. Michael On 10/01/2023 12:30, Fahimeh Alizadeh via R-help wrote: Dear Prof. My forest plot is not fit to windows in R software, I have 49 studies but forest plot only show from 9 to 42. forest.meta(hidemeta, layout="RevMan5", xlab="Proportion", comb.r=T, comb.f=F, xlim = c(0,1), fontsize=10, digits=3) Whould you please help me to access full forest plot? Thank you Sincerely,Fahimeh Dr. Fahimeh AlizadehPost-Doctorate Fellow and Ph.DMicrobial Biotechnology Islamic Azad University __ 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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] Functional Programming Problem Using purr and R's data.table shift function
Dénes, thank you for the guidance - which is well-taken.
Your side note raises an interesting question: I find the piping %>%
operator readable. Is there any downside to it? Or is the side note meant
to tell me to drop the last: "%>% `[`"?
Thank you,
==
Michael Lachanski
PhD Student in Demography and Sociology
MA Candidate in Statistics
University of Pennsylvania
mikel...@sas.upenn.edu
On Sat, Dec 31, 2022 at 9:22 AM Dénes Tóth wrote:
> Hi Michael,
>
> Note that you have to be very careful when using by-reference operations
> in data.table (see `?data.table::set`), especially in a functional
> programming approach. In your function, you avoid this problem by
> calling `data.table(A)` which makes a copy of A even if it is already a
> data.table. However, for large data.table-s, copying can be a very
> expensive operation (esp. in terms of RAM usage), which can be totally
> eliminated by using data.tables in the data.table-way (e.g., joining,
> grouping, and aggregating in the same step by performing these
> operations within `[`, see `?data.table`).
>
> So instead of blindly functionalizing all your code, try to be
> pragmatic. Functional programming is not about using pure functions in
> *every* part of your code base, because it is unfeasible in 99.9% of
> real-world problems. Even Haskell has `IO` and `do`; the point is that
> the imperative and functional parts of the code are clearly separated
> and imperative components are (tried to be) as top-level as possible.
>
> So when using data.table, a good strategy is to use pure functions for
> performing within-data.table operations, e.g., `DT[, lapply(.SD, mean),
> .SDcols = is.numeric]`, and when these operations alter `DT` by
> reference, invoke the chains of these operations in "pure" wrappers -
> e.g., calling `A <- copy(A)` on the top and then modifying `A` directly.
>
> Cheers,
> Denes
>
> Side note: You do not need to use `DT[ , A:= shift(A, fill = NA, type =
> "lag", n = 1)] %>% `[`(return(DT))`. `[.data.table` returns the result
> (the modified DT) invisibly. If you want to let auto-print work, you can
> just use `DT[ , A:= shift(A, fill = NA, type = "lag", n = 1)][]`.
>
> Note that this also means you usually you do not need to use magrittr's
> or base-R pipe when transforming data.table-s. You can do this instead:
> ```
> DT[
>## filter rows where 'x' column equals "a"
>x == "a"
> ][
>## calculate the mean of `z` for each gender and assign it to `y`
>, y := mean(z), by = "gender"
> ][
>## do whatever you want
>...
> ]
> ```
>
>
> On 12/31/22 13:39, Rui Barradas wrote:
> > Às 06:50 de 31/12/2022, Michael Lachanski escreveu:
> >> Hello,
> >>
> >> I am trying to make a habit of "functionalizing" all of my code as
> >> recommended by Hadley Wickham. I have found it surprisingly difficult
> >> to do
> >> so because several intermediate features from data.table break or give
> >> unexpected results using purrr and its data.table adaptation, tidytable.
> >> Here is the a minimal working example of what has stumped me most
> >> recently:
> >>
> >> ===
> >>
> >> library(data.table); library(tidytable)
> >>
> >> minimal_failing_function <- function(A){
> >>DT <- data.table(A)
> >>DT[ , A:= shift(A, fill = NA, type = "lag", n = 1)] %>% `[`
> >>return(DT)}
> >> # works
> >> minimal_failing_function(c(1,2))
> >> # fails
> >> tidytable::pmap_dfr(.l = list(c(1,2)),
> >> .f = minimal_failing_function)
> >>
> >>
> >> ===
> >> These should ideally give the same output, but do not. This also fails
> >> using purrr::pmap_dfr rather than tidytable. I am using R 4.2.2 and I
> >> am on
> >> Mac OS Ventura 13.1.
> >>
> >> Thank you for any help you can provide or general guidance.
> >>
> >>
> >> ==
> >> Michael Lachanski
> >> PhD Student in Demography and Sociology
> >> MA Candidate in Statistics
> >> University of Pennsylvania
> >> mikel...@sas.upenn.edu
> >>
> >> [[alternative HTML version deleted]]
> >>
> >> __
> >> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> >>
> https://urldefense.com/v3/__https://stat.ethz.ch/mailman/listinfo/r-help__;!!IBzWLUs!VdfzdJ15GLScUok_hiqL3DvTJ20Ce8JMBkQ1NosBfyOvu68iuQkh9nsPZuUBbB9BtrsZBh86OjGyyj3lAB2g_xXCvB6
[R] Functional Programming Problem Using purr and R's data.table shift function
Hello,
I am trying to make a habit of "functionalizing" all of my code as
recommended by Hadley Wickham. I have found it surprisingly difficult to do
so because several intermediate features from data.table break or give
unexpected results using purrr and its data.table adaptation, tidytable.
Here is the a minimal working example of what has stumped me most recently:
===
library(data.table); library(tidytable)
minimal_failing_function <- function(A){
DT <- data.table(A)
DT[ , A:= shift(A, fill = NA, type = "lag", n = 1)] %>% `[`
return(DT)}
# works
minimal_failing_function(c(1,2))
# fails
tidytable::pmap_dfr(.l = list(c(1,2)),
.f = minimal_failing_function)
===
These should ideally give the same output, but do not. This also fails
using purrr::pmap_dfr rather than tidytable. I am using R 4.2.2 and I am on
Mac OS Ventura 13.1.
Thank you for any help you can provide or general guidance.
==
Michael Lachanski
PhD Student in Demography and Sociology
MA Candidate in Statistics
University of Pennsylvania
mikel...@sas.upenn.edu
[[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] piecewiseSEM two errors: arguments imply different numbers of rows / object 'ret' not found
Dear Rlist members,
We are trying to run a SEM with your package. Unfortunaley while running our
code two error messages popped up. For the life of us we could not figure
out how to solve these errors . We are hoping very much that somebody can
help us: Why do we get these error messages? We built our model structure
using an example from the piecewiseSEM package and we double checked the
that we do not have any "NAs" in our data set. Thanks a ton!
Error1:
Error in data.frame(..., check.names = FALSE) :
arguments imply different numbers of rows: 166, 0
Error2: (we ran the model with a smaller,simpler SEM)
Error in cbind(ret, isSig(ret[, 5])) : object 'ret' not found
Below is our code and the data (as "dput()")
# dataCode: (see dput(dataCode) at the very end
library (glmmTMB)
library(lme4)
library(DHARMa)#check for Overdispersion
library(piecewiseSEM)
library(lavaan)
str(dataCode)
summary(dataCode)
#define factors
dataCode$Year<-factor(dataCode$Year)
#relevel(dataCode$Year,ref="2019")
dataCode$Stand<-as.factor(dataCode$Stand)
dataCode$TreeNo <-as.factor(dataCode$TreeNo)
dataCode$Drought <-as.factor(dataCode$Drought)
#relevel(dataCode$Drought,ref="0")
dataCode$Stratum <-as.factor(dataCode$Stratum)
levels(dataCode$Stratum) <- list(lower = "shade", upper = "sun")
#relevel(dataCode$Stratum,ref="sun")
dataCode$Location <-as.factor(dataCode$Location)
dataCode$Tree.ID <-as.factor(dataCode$Tree.ID)
# Models for SEM
mod1x<-lmer(N_pc~StratumDrought+YearDrought+(1|Tree.ID),data)
#Check assumptions
simulationOutput <- simulateResiduals(fittedModel = mod1x)
plot(simulationOutput)
summary(mod1x)
mod2x<-lmer(log(Fiber)~StratumDrought+YearDrought+(1|Tree.ID),data)
#Check assumptions
simulationOutput <- simulateResiduals(fittedModel = mod2x)
plot(simulationOutput)
summary(mod2x)
mod3x<-lmer(Lignin~StratumDrought+YearDrought+(1|Tree.ID),data)
#Check assumptions
simulationOutput <- simulateResiduals(fittedModel = mod3x)
plot(simulationOutput)
summary(mod3x)
mod11a<-glmer(prop_suck~N_pc+(1|Tree.ID),
data=data, family=binomial(link = "logit"), weights=weight_suck)
mod12a<-glmmTMB(prop_suck~Fiber+(1|Tree.ID),
data=data, family=betabinomial(link = "logit"), weights=weight_suck)
mod13a<-glmmTMB(prop_suck~Lignin+(1|Tree.ID),
data=data, family=betabinomial(link = "logit"), weights=weight_suck)
newlist2 = list(
mod1x, mod2x, mod3x, mod11a, mod12a, mod13a)
model<-as.psem(newlist2)#Error1:---!
#We then ran a model with a simpler model structure, to double check if the
error occured due to the model structure
#Error 1 did not appear anymore. However, error 2 popped up-!
modlist = list(
mod1x, mod11a)
model<-as.psem(modlist)
summary(model) #Error1:
#DATA
dput(dataCode)
structure(list(Tree.ID = structure(c(1L, 1L, 2L, 2L, 3L, 3L,
4L, 4L, 5L, 5L, 6L, 6L, 7L, 7L, 8L, 8L, 9L, 9L, 10L, 10L, 11L,
11L, 12L, 12L, 13L, 13L, 14L, 14L, 15L, 15L, 16L, 16L, 17L, 17L,
18L, 18L, 19L, 19L, 20L, 20L, 21L, 21L, 22L, 22L, 23L, 23L, 24L,
24L, 25L, 25L, 26L, 26L, 27L, 27L, 28L, 28L, 29L, 29L, 30L, 30L,
31L, 31L, 32L, 32L, 33L, 33L, 34L, 34L, 35L, 35L, 36L, 36L, 37L,
37L, 38L, 38L, 39L, 39L, 40L, 40L, 41L, 41L, 42L, 42L, 43L, 43L,
44L, 44L, 45L, 45L, 46L, 46L, 47L, 47L, 48L, 48L, 1L, 1L, 2L,
2L, 3L, 3L, 10L, 10L, 11L, 11L, 12L, 12L, 13L, 13L, 14L, 14L,
15L, 15L, 16L, 16L, 17L, 17L, 18L, 18L, 19L, 19L, 21L, 21L, 22L,
22L, 23L, 23L, 24L, 24L, 25L, 25L, 26L, 26L, 27L, 27L, 28L, 28L,
29L, 29L, 30L, 30L, 31L, 31L, 32L, 32L, 33L, 33L, 37L, 37L, 38L,
38L, 39L, 39L, 43L, 43L, 44L, 44L, 45L, 45L, 46L, 46L, 47L, 47L,
48L, 48L), .Label = c("102_6", "102_7", "102_8", "105_1", "105_2",
"105_4", "111_7", "111_8", "111_9", "113_2", "113_4", "113_5",
"114_7", "114_8", "114_9", "116_6", "116_7", "116_9", "122_3",
"122_4", "122_5", "132_3", "132_4", "132_5", "242_2", "242_4",
"242_5", "243_1", "243_2", "243_4", "245_1", "245_2", "245_5",
"246_1", "246_2", "246_3", "251_10", "251_8", "251_9", "253_7",
"253_8", "253_9", "254_6", "254_7", "254_8", "267_10", "267_6",
"267_8"), class = "factor"), Year = structure(c(1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L), .Label = c("1", "2"), class = "factor"), Drought =
structure(c(1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
Re: [R] Select dataframe row containing a digit
I suspect [[:digit:]] might have done what you want.
Michael
On 30/11/2022 13:04, Luigi Marongiu wrote:
Thank you,
I have been trying with [:digit:] but did not work. It worked with
`df$val[grepl('[0-9]', df$val)] = "NUM"`
On Wed, Nov 30, 2022 at 2:02 PM Ivan Krylov wrote:
В Wed, 30 Nov 2022 13:40:50 +0100
Luigi Marongiu пишет:
I am formatting everything to either "POS" and "NEG",
but values entered as number should get the value "NUM".
How do I change such values?
Thanks for providing an example!
One idea would be to use a regular expression to locate numbers. For
example, grepl('[0-9]', df$val) will return a logical vector indexing
the rows containing digits. Alternatively, grepl('^[0-9.]+$', df$val,
perl = TRUE) will index all strings consisting solely of digits and
decimal separators.
Another idea would be to parse all of the strings as numbers and filter
out those that didn't succeed. Use as.numeric() to perform the parsing,
suppressWarnings() to silence the messages telling you that the parsing
failed for some of the strings and is.na() to get the logical vector
indexing those entries that failed to parse.
--
Best regards,
Ivan
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] picewiseSEM undefined columns selected
Dear r-help list members,
I would like to run a pathway analysis using the R-package piecewiseSEM
(v.2.12). I used a very simple model structure that is similar to the example
in the piecewiseSEM description. Unfortunatley, once I execute the "as.psem"
function I get the error "Error in `[.data.frame`(x$data, , vars) : undefined
columns selected" . This souns like a simple error message. Yet, I checked my
data and the model structure several times and I can't figure out why the error
appears. Does anybody know why I still get the error? I posted a simple code
example and the data set that goes with it below.
Many thanks in advance
Michi
CODE:
data2$Year<-as.factor(data2$Year)
data2$Drought<-as.factor(data2$Drought)
data2$Stratum<-as.factor(data2$Stratum)
mod1x<-lmer(N_pc~Stratum:Drought+Year:Drought+(1|Tree.ID),data2)
mod11a<-glmer(cbind(branch_suck_Yes,branch_suck_No)~N_pc+(1|Tree.ID),
data=data2, family = binomial(link = "logit"))
modlist = list(
mod1x, mod11a)
model<-as.psem(modlist)#Error: Error in `[.data.frame`(x$data, , vars) :
undefined columns selected
model<-psem(modlist, data2)
coefs(modlist, data2, standardize = "none", intercept = FALSE)
model<-psem(mod1x, mod11a)
DATA:
structure(list(N_pc = c(2.39, 2.81, 2.48, 1.83, 1.91, 1.96, 2.32,
2.54, 2.19, 1.97, 2.29, 1.64, 2.03, 1.68, 2.145, 2.08, 1.99,
1.99, 2.83, 2.83, 2.91, 2.61, 2.73, 2.54, 2.87, 1.91, 2.84, 2.74,
2.87, 2.6, 2.12, 2.64, 2.46, 1.83, 2.06, 2.77, 2.41, 2.74, 2.83,
2.51, 2.79, 2.66, 2.44, 2.26, 2.85, 2.39, 2.52, 2.13, 2.63, 2,
2.43, 2.36, 2.98, 2.28, 2.12, 2.2, 2.54, 1.28, 2.57, 2.17, 2.32,
2.41, 3.11, 2.591, 2.77, 2.53, 2.67, 2.45, 2.5, 2.52, 2.9, 3.03,
2.83, 2.52, 2.57, 2.62, 2.82, 2.62, 2.98, 3.01, 2.33, 2.11, 2.68,
2.74, 2.53, 2.43), Stratum = structure(c(1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 2L,
1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L), .Label = c("shade", "sun"), class = "factor"), Drought = structure(c(1L,
2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L), .Label = c("no", "yes"), class = "factor"),
Year = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L), .Label = c("1", "2"), class = "factor"), Tree.ID = structure(c(2L,
15L, 19L, 21L, 22L, 23L, 25L, 28L, 29L, 29L, 30L, 31L, 32L,
33L, 34L, 34L, 35L, 36L, 1L, 1L, 2L, 2L, 3L, 3L, 4L, 4L,
5L, 5L, 6L, 6L, 7L, 7L, 8L, 8L, 9L, 10L, 10L, 11L, 11L, 12L,
13L, 13L, 14L, 14L, 15L, 15L, 16L, 16L, 17L, 17L, 18L, 18L,
19L, 19L, 20L, 20L, 21L, 21L, 22L, 22L, 23L, 23L, 24L, 24L,
25L, 25L, 26L, 26L, 27L, 27L, 28L, 28L, 29L, 29L, 32L, 32L,
33L, 33L, 34L, 34L, 35L, 35L, 36L, 36L, 37L, 37L), .Label = c("102_6",
"102_7", "102_8", "113_2", "113_4", "113_5", "114_7", "114_8",
"114_9", "116_6", "116_7", "116_9", "122_3", "122_5", "132_3",
"132_4", "132_5", "242_2", "242_4", "242_5", "243_1", "243_2",
"243_4", "245_1", "245_2", "245_5", "251_10", "251_8", "251_9",
"253_7", "253_8", "254_6", "254_7", "254_8", "267_10", "267_6",
"267_8"), class = "factor"), branch_suck_Yes = c(2L, 3L,
1L, 6L, 6L, 5L, 4L, 8L, 2L, 2L, 2L, 48L, 3L, 22L, 14L, 2L,
1L, 1L, 27L, 25L, 16L, 13L, 31L, 18L, 31L, 2L, 25L, 16L,
36L, 21L, 8L, 23L, 13L, 11L, 7L, 35L, 7L, 17L, 4L, 40L, 48L,
17L, 16L, 10L, 34L, 6L, 46L, 7L, 26L, 12L, 24L, 27L, 31L,
25L, 34L, 21L, 44L, 23L, 40L, 30L, 25L, 18L, 35L, 10L, 6L,
10L, 29L, 5L, 24L, 16L, 19L, 11L, 21L, 10L, 18L, 18L, 42L,
33L, 16L, 31L, 16L, 38L, 37L, 28L, 9L, 35L), branch_suck_No = c(48L,
47L, 49L, 44L, 44L, 45L, 46L, 42L, 48L, 48L, 48L, 2L, 47L,
28L, 36L, 48L, 49L, 49L, 23L, 25L, 34L, 37L, 19L, 32L, 19L,
48L, 25L, 34L, 14L, 29L, 42L, 27L, 37L, 39L, 43L, 15L, 43L,
33L, 46L, 10L, 2L, 33L, 34L, 40L, 16L, 44L, 4L, 43L, 24L,
38L, 26L, 23L, 19L, 25L, 16L, 29L, 6L, 27L, 10L, 20L, 25L,
32L, 15L, 40L, 44L, 40L, 21L, 45L, 26L, 34L, 31L, 39L, 29L,
40L, 32L, 32L, 8L, 17L, 34L, 19L, 34L, 12L, 13L, 22L, 41L,
15L)), row.names = c(3L, 43L, 51L, 55L, 57L, 59L, 63L, 76L,
77L, 78L, 80L, 81L, 86L, 87L, 89L, 90L, 92L, 93L, 97L, 98L, 99L,
100L, 101L, 102L, 103L, 104L, 105L, 106L, 107L, 108L, 109L, 110L,
111L, 112L, 113L, 115L, 116L, 117L, 118L, 119L, 121L, 122L, 123L,
124L, 125L, 126L, 127L, 128L, 129L, 130L, 131L, 132L, 133L, 134L,
135L,
Re: [R] Sensitivity and Specificity
So predict is a one-dimensional vector of predictions but you are
treating it as a two-dimensional matrix (presumably you think those are
the totals).
Michael
On 24/10/2022 16:50, greg holly wrote:
Hi Michael,
I appreciate your writing. Here are what I have after;
> predict_testing <- ifelse(predict > 0.5,1,0)
>
> head(predict)
1 2 3 5 7 8
0.29006984 0.28370507 0.10761993 0.02204224 0.12873872 0.08127920
>
> # Sensitivity and Specificity
>
>
sensitivity<-(predict_testing[2,2]/(predict_testing[2,2]+predict_testing[2,1]))*100
Error in predict_testing[2, 2] : incorrect number of dimensions
> sensitivity
function (data, ...)
{
UseMethod("sensitivity")
}
>
>
specificity<-(predict_testing[1,1]/(predict_testing[1,1]+predict_testing[1,2]))*100
Error in predict_testing[1, 1] : incorrect number of dimensions
> specificity
function (data, ...)
{
UseMethod("specificity")
}
On Mon, Oct 24, 2022 at 10:45 AM Michael Dewey <mailto:li...@dewey.myzen.co.uk>> wrote:
Rather hard to know without seeing what output you expected and what
error message you got if any but did you mean to summarise your
variable
predict before doing anything with it?
Michael
On 24/10/2022 16:17, greg holly wrote:
> Hi all R-Help ,
>
> After partitioning my data to testing and training (please see
below), I
> need to estimate the Sensitivity and Specificity. I failed. It
would be
> appropriate to get your help.
>
> Best regards,
> Greg
>
>
> inTrain <- createDataPartition(y=data$case,
> p=0.7,
> list=FALSE)
> training <- data[ inTrain,]
> testing <- data[-inTrain,]
>
> attach(training)
> #model training and prediction
> data_training <- glm(case ~ age+BMI+Calcium+Albumin+meno_1, data =
> training, family = binomial(link="logit"))
>
> predict <- predict(data_training, data_predict = testing, type =
"response")
>
> predict_testing <- ifelse(predict > 0.5,1,0)
>
> # Sensitivity and Specificity
>
>
sensitivity<-(predict_testing[2,2]/(predict_testing[2,2]+predict_testing[2,1]))*100
> sensitivity
>
>
specificity<-(predict_testing[1,1]/(predict_testing[1,1]+predict_testing[1,2]))*100
> specificity
>
> [[alternative HTML version deleted]]
>
> __
> R-help@r-project.org <mailto:R-help@r-project.org> mailing list
-- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
<https://stat.ethz.ch/mailman/listinfo/r-help>
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
<http://www.R-project.org/posting-guide.html>
> and provide commented, minimal, self-contained, reproducible code.
>
--
Michael
http://www.dewey.myzen.co.uk/home.html
<http://www.dewey.myzen.co.uk/home.html>
<http://www.avg.com/email-signature?utm_medium=email_source=link_campaign=sig-email_content=emailclient>
Virus-free.www.avg.com
<http://www.avg.com/email-signature?utm_medium=email_source=link_campaign=sig-email_content=emailclient>
<#DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2>
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] Sensitivity and Specificity
Rather hard to know without seeing what output you expected and what error message you got if any but did you mean to summarise your variable predict before doing anything with it? Michael On 24/10/2022 16:17, greg holly wrote: Hi all R-Help , After partitioning my data to testing and training (please see below), I need to estimate the Sensitivity and Specificity. I failed. It would be appropriate to get your help. Best regards, Greg inTrain <- createDataPartition(y=data$case, p=0.7, list=FALSE) training <- data[ inTrain,] testing <- data[-inTrain,] attach(training) #model training and prediction data_training <- glm(case ~ age+BMI+Calcium+Albumin+meno_1, data = training, family = binomial(link="logit")) predict <- predict(data_training, data_predict = testing, type = "response") predict_testing <- ifelse(predict > 0.5,1,0) # Sensitivity and Specificity sensitivity<-(predict_testing[2,2]/(predict_testing[2,2]+predict_testing[2,1]))*100 sensitivity specificity<-(predict_testing[1,1]/(predict_testing[1,1]+predict_testing[1,2]))*100 specificity [[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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] Forestplot, grid graphics Viewplot grid.arange
There are several other CRAN packages which provide forest plots (see
CRAN Task View for details) and they do not all use grip graphics which
I think the forestplot package does. It might be worth swapping to one
of them.
Michael
On 14/10/2022 04:34, Jim Lemon wrote:
Hi Mary,
I didn't see any answers to your post, but doing something like this
is quite easy in base graphics. If you are still stuck, I may be able
to suggest something.
Jim
On Mon, Oct 10, 2022 at 6:05 PM Putt, Mary wrote:
I have created several plots using the forestplot package and the link shown here.
<https://cran.r-project.org/web/packages/forestplot/vignettes/forestplot.html >
Great package !
Next step is to combine two plots into a single graphic. The code provided on
the link results in 'bleeding' of the graphics/text into each other. I don't
want to clip it as I need the text elements. I am guessing the problem involves
the combination of text and graphics in the 'plot'. I fooled around with the
original post and also did some hunting online but no luck Thanks in advance.
library(foresplot)
data("dfHRQoL")
#create individual forest plots for Sweden and Denmark
fp_sweden <- dfHRQoL |>
filter(group == "Sweden") |>
mutate(est = sprintf("%.2f", mean), .after = labeltext) |>
forestplot(labeltext = c(labeltext, est),
title = "Sweden",
clip = c(-.1, Inf),
xlab = "EQ-5D index",
new_page = FALSE)
fp_denmark <- dfHRQoL |>
filter(group == "Denmark") |>
mutate(est = sprintf("%.2f", mean), .after = labeltext) |>
forestplot(labeltext = c(labeltext, est),
title = "Denmark",
clip = c(-.1, Inf),
xlab = "EQ-5D index",
new_page = FALSE)
#now combine into a single plot using the web code; but this one bleeds into
each other
library(grid)
#
#Put plots together using grid graphics
#Attempt 1 from website
#
grid.newpage()
borderWidth <- unit(4, "pt")
width <- unit(convertX(unit(1, "npc") - borderWidth, unitTo = "npc", valueOnly = TRUE)/2,
"npc")
pushViewport(viewport(layout = grid.layout(nrow = 1,
ncol = 3,
widths = unit.c(width,
borderWidth,
width))
)
)
pushViewport(viewport(layout.pos.row = 1,
layout.pos.col = 1))
fp_sweden
upViewport()
pushViewport(viewport(layout.pos.row = 1,
layout.pos.col = 2))
grid.rect(gp = gpar(fill = "grey", col = "red"))
upViewport()
pushViewport(viewport(layout.pos.row = 1,
layout.pos.col = 3))
fp_denmark
upViewport(2)
#Attempt 2 from website, still a problem.
grid.newpage()
borderWidth <- unit(4, "pt")
width <- unit(convertX(unit(1, "npc") - borderWidth, unitTo = "npc", valueOnly = TRUE)/2,
"npc")
pushViewport(viewport(layout = grid.layout(nrow = 1,
ncol = 3,
widths = c(0.45, 0.1, 0.45))
)
)
pushViewport(viewport(layout.pos.row = 1,
layout.pos.col = 1))
fp_sweden
upViewport()
pushViewport(viewport(layout.pos.row = 1,
layout.pos.col = 3))
fp_denmark
upViewport(2)
###
#Attempt 3 converting to grobs and use patchwork
###
library(ggplotify)
library(patchwork)
fpd_grob <- grid2grob(print(fp_denmark))
p1 <- grid2grob(print(fp_denmark))
p2 <- grid2grob(print(fp_sweden))
p_both <- wrap_elements(p1) + wrap_elements(p2)
p_both
#same problem with grid.arrange()**strong text**
Mary Putt, PhD, ScD
Professor of Biostatistics
Department of Biostatistics, Epidemiology & Informatics
Pereleman School of Medicine
University of Pennsylvania
215-573-7020
[[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.
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] High concurvity/ collinearity between time and temperature in GAM predicting deaths but low ACF. Does this matter?
Dear Jade
It seems to me that there are several issues here.
1 - if you fit a separate parameter for each heaping day you efectively
remove them from the model altogether. If they do carry meaningful
information then that is undesirable.
2 - if the reason that there ae so many 15s is because people state,
"Oh, it was in July" and either they or the interviewere impute 15
because it si the month mid-point then it begins to look as though
attempting to use daily data is fraught with problems and monthly might
be better.
3 - if you fit a binary heap variable which is zero except for the 15 of
the month and 1 for the 15 then you are taking out any general overall
tendency to report 15 but you are allowing for the possibility that the
relative size of heaping days can vary during the study period.
I am not an expert on this sort of model, that is Simon clearly, but it
seems to me that, as usual, we have a mix of statistical and scientific
questions here and you may need to rethink your sceientific model in the
light of the issue with your data.
Michael
On 10/06/2022 10:30, jade.sho...@googlemail.com wrote:
Hi Michael,
I don't think my reply to your email came through to the list, so am
resending (see below). Problems with subscription have now hopefully
been resolved. Apologies if this is a double posting!
On Thu, 9 Jun 2022 at 15:27, jade.sho...@googlemail.com
wrote:
Hi Michael,
Thanks for the reply! When I ran the gam with the gam() function, the
model worked fine with heap having 169 levels. The same model with
bam() however, fails. I don't understand the difference between bam()
and gam() at all (other than computational efficiency), but could the
fact that each level only has 1 data point be the reason for it?
These heaping days are very large measurement errors I need to get rid
off, so I just want to take them out of the model altogether. (My data
is quite noisy already, because date of death is based on memory
recall, rather than formal death registration, due to the data being
from a low income country). Median of deaths is approx. 2 per day, but
on heaping days it can be as high as 50 or so.
My understanding was that coding with 169 levels would effectively
take these measurements out of the model (but do correct me if I'm
wrong!) I originally coded heap as a binary variable with 0 for
non-heaping days and 1 for heaping days, but was told that meant that
I was assuming the effect was the same for all heaping days. If I
coded with 12 or 14 levels, wouldn't that leave a lot of noise in the
data?
Jade
On Thu, 9 Jun 2022 at 14:52, Michael Dewey wrote:
Dear Jade
Do you really need to fit a separate parameter for each heaping day? Can
you not just make it a binary predictor or a categorical one with fewer
levels, perhaps 14 (for heaping in each year) or 12 (for each calendar
month). I have no idea whether that would help but it seems worth a try.
Michael
On 08/06/2022 18:15, jade.shodan--- via R-help wrote:
Hi Simon,
Thanks so much for this!! I have two follow up questions, if you don't mind.
1. Does including an autoregressive term not adjust away part of the
effect of the response in a distributed lag model (where the outcome
accumulates over time)?
2. I've tried to fit a model using bam (just a first attempt without
AR term), but including the factor variable heap creates errors:
bam0 <- bam(deaths~te(year, month, week, weekday,
bs=c("cr","cc","cc","cc"), k = c(4,5,5,5)) + heap +
te(temp_max, lag, k=c(8, 3)) +
te(precip_daily_total, lag, k=c(8, 3)),
data = dat, family = nb, method = 'fREML',
select = TRUE, discrete = TRUE,
knots = list(month = c(0.5, 12.5), week = c(0.5,
52.5), weekday = c(0, 6.5)))
This model results in errors:
Warning in estimate.theta(theta, family, y, mu, scale = scale1, wt = G$w, :
step failure in theta estimation
Warning in sqrt(family$dev.resids(object$y, object$fitted.values,
object$prior.weights)) :
NaNs produced
Including heap as as.numeric(heap) runs the model without error
messages or warnings, but model diagnostics look terrible, and it also
doesn't make sense (to me) to make heap a numeric. The factor variable
heap (with 169 levels) codes the fact that all deaths for which no
date was known, were registered on the 15th day of each month. I've
coded all non-heaping days as 0. All heaping days were coded as a
value between 1-168. The time series spans 14 years, so a heaping day
in each month results in 14*12 levels = 168, plus one level for
non-heaping days.
So my second question is: Does bam allow factor variables? And if not,
how should I model this heaping on the 15th day of the month instead?
With thanks,
Jade
On Wed, 8 Jun 2022 at 12:05, Simon Wood wrote:
I would not worry too much about high concurvity between variables like
temperature and time. This
Re: [R] High concurvity/ collinearity between time and temperature in GAM predicting deaths but low ACF. Does this matter?
Dear Jade
Do you really need to fit a separate parameter for each heaping day? Can
you not just make it a binary predictor or a categorical one with fewer
levels, perhaps 14 (for heaping in each year) or 12 (for each calendar
month). I have no idea whether that would help but it seems worth a try.
Michael
On 08/06/2022 18:15, jade.shodan--- via R-help wrote:
Hi Simon,
Thanks so much for this!! I have two follow up questions, if you don't mind.
1. Does including an autoregressive term not adjust away part of the
effect of the response in a distributed lag model (where the outcome
accumulates over time)?
2. I've tried to fit a model using bam (just a first attempt without
AR term), but including the factor variable heap creates errors:
bam0 <- bam(deaths~te(year, month, week, weekday,
bs=c("cr","cc","cc","cc"), k = c(4,5,5,5)) + heap +
te(temp_max, lag, k=c(8, 3)) +
te(precip_daily_total, lag, k=c(8, 3)),
data = dat, family = nb, method = 'fREML',
select = TRUE, discrete = TRUE,
knots = list(month = c(0.5, 12.5), week = c(0.5,
52.5), weekday = c(0, 6.5)))
This model results in errors:
Warning in estimate.theta(theta, family, y, mu, scale = scale1, wt = G$w, :
step failure in theta estimation
Warning in sqrt(family$dev.resids(object$y, object$fitted.values,
object$prior.weights)) :
NaNs produced
Including heap as as.numeric(heap) runs the model without error
messages or warnings, but model diagnostics look terrible, and it also
doesn't make sense (to me) to make heap a numeric. The factor variable
heap (with 169 levels) codes the fact that all deaths for which no
date was known, were registered on the 15th day of each month. I've
coded all non-heaping days as 0. All heaping days were coded as a
value between 1-168. The time series spans 14 years, so a heaping day
in each month results in 14*12 levels = 168, plus one level for
non-heaping days.
So my second question is: Does bam allow factor variables? And if not,
how should I model this heaping on the 15th day of the month instead?
With thanks,
Jade
On Wed, 8 Jun 2022 at 12:05, Simon Wood wrote:
I would not worry too much about high concurvity between variables like
temperature and time. This just reflects the fact that temperature has a
strong temporal pattern.
I would also not be too worried about the low p-values on the k check.
The check only looks for pattern in the residuals when they are ordered
with respect to the variables of a smooth. When you have time series
data and some smooths involve time then it's hard not to pick up some
degree of residual auto-correlation, which you often would not want to
model with a higher rank smoother.
The NAs for the distributed lag terms just reflect the fact that there
is no obvious way to order the residuals w.r.t. the covariates for such
terms, so the simple check for residual pattern is not really possible.
One simple approach is to fit using bam(...,discrete=TRUE) which will
let you specify an AR1 parameter to mop up some of the residual
auto-correlation without resorting to a high rank smooth that then does
all the work of the covariates as well. The AR1 parameter can be set by
looking at the ACF of the residuals of the model without this. You need
to look at the ACF of suitably standardized residuals to check how well
this has worked.
best,
Simon
On 05/06/2022 20:01, jade.shodan--- via R-help wrote:
Hello everyone,
A few days ago I asked a question about concurvity in a GAM (the
anologue of collinearity in a GLM) implemented in mgcv. I think my
question was a bit unfocussed, so I am retrying again, but with
additional information included about the autocorrelation function. I
have also posted about this on Cross Validated. Given all the model
output, it might make for easier
reading:https://stats.stackexchange.com/questions/577790/high-concurvity-collinearity-between-time-and-temperature-in-gam-predicting-dea
As mentioned previously, I have problems with concurvity in my thesis
research, and don't have access to a statistician who works with time
series, GAMs or R. I'd be very grateful for any (partial) answer,
however short. I'll gladly return the favour where I can! For really
helpful input I'd be more than happy to offer co-authorship on
publication. Deadlines are very close, and I'm heading towards having
no results at all if I can't solve this concurvity issue :(
I'm using GAMs to try to understand the relationship between deaths
and heat-related variables (e.g. temperature and humidity), using
daily time series over a 14-year period from a tropical, low-income
country. My aim is to understand the relationship between these
variables and deaths, rather than pure prediction performance.
The GAMs include distributed lag models (set up as 7-column matrices,
see code at bottom of post), since deaths may occur over several days
foll
Re: [R] Shiny app Deployment not Working
I think people are going to need more details here. 1 where does GeneBook come from if not CRAN? 2 what is the app you are hoping to deploy, is it part of GeneBook? 3 what exactly do you mean by not showing? On 27/05/2022 07:29, ABHIRUP MOITRA wrote: I want to deploy a shiny app where I have to use the package 'GeneBook' which is not available in CRAN. So the app is not showing after deployment. Thanks and regards Abhirup Moitra [[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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] Restoration of "rite" package of R as R-script editor or the like
Notepad++ is fine as a light weight editor with syntax highlighting but the plug-in for R does not seem to have been updated for years. Michael On 23/05/2022 17:37, Rui Barradas wrote: Hello, I was also thinking about Notepad++. And about an editor no one mentinoed, vim. I've learned vi on unix SVR3 and at least in micros it is all present, you can bet that no matter what the OS is there's a version for it. vim is lightweight and has syntax highlighting but its aging idiosyncrasies might outweight its benefits. Rui Barradas Às 07:35 de 22/05/2022, Ivan Krylov escreveu: Dear Dr. A.K. Singh, The person responsible for the "rite" package is Dr. Thomas J. Leeper, and he doesn't seem interested in maintaining it further: https://github.com/leeper/rite (I do like the idea of a Tcl/Tk text editor with tight R integration, but maintaining it, like maintaining anything, is a job one needs money and/or other motivation to perform.) Have you considered using some other programmer's text editor with R, like Notepad++? I've heard there's a companion utility that helps integrate R and Notepad++: https://github.com/halpo/NppToR __ 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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] How to create density ellipses with R
ellipse draws mathematical ellipses, I believe For proper bivariate normal density ellipses (not simply contours), check out: car::dataEllipse heplots::covEllipses -Michael Message: 1 Date: Fri, 14 Jan 2022 10:12:50 -0500 From: Paul Bernal To: R Subject: [R] How to create density ellipses with R Message-ID: Content-Type: text/plain; charset="utf-8" Dear R friends, Happy new year to you all. Not quite sure if this is the proper place to ask about this, so I apologize if it is not, and if it isn´t, maybe you can point me to the right place. I would like to know if there is any R package that allows me to produce density ellipses. Searching through the net, I came across a package called ellipse, but I'm not sure if this is the one I should use. Any help and/or guidance will be greatly appreciated. Best regards, Paul [[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: [ESS] ESS process on Docker containers and enabing Flymake
To be clear, Alex is referring to the R package 'lintr', not an emacs package. I made that mistake and spent some time looking for an elpa package. Mike __ ESS-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/ess-help
Re: [R] how to run r biotools boxM terst on multiple groups?
Try heplots::boxM This also has a nice plot method to visualize the contributions to box's M test. See my paper: https://www.datavis.ca/papers/EqCov-TAS.pdf > -Original Message- > From: R-help > mailto:r-help-boun...@r-project.org>> On Behalf > Of Luigi > Marongiu > Sent: Tuesday, January 4, 2022 11:56 AM > To: r-help mailto:r-help@r-project.org>> > Subject: [R] how to run r biotools boxM terst on multiple groups? > > I have a data frame containing a half dozen continuous measurements > and over a dozen ordinal variables (such as, death, fever, symptoms etc). > I would like to run a box matrix test and I am using biotools' boxM, > but it Michael Friendly Email: friendly AT yorku DOT ca Professor, Psychology Dept. & Former Chair, ASA Statistical Graphics Section York University Voice: 416 736-2100 x66249 4700 Keele StreetWeb: http://www.datavis.ca | @datavisFriendly Toronto, ONT M3J 1P3 CANADA [[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] Issue about eRM package
I do not know that package so cannot help with that but the dataset did not come through. This mailing list is quite restrictive in what sorts of attachment it allows so I suggest trying something like a .csv file or a .txt file. Michael On 28/09/2021 10:37, 宋启发 wrote: I am using your excellent R package eRM to solve a questionaire survey data. I meet a strange issue when using RSM function. This RSM runs well using sample data. There is a test data like this: 'data.frame': 1669 obs. of 7 variables: $ X1: num 2 2 3 3 2 3 4 4 3 2 ... $ X2: num 4 3 3 2 3 4 4 4 3 3 ... $ X3: num 4 3 3 4 3 3 4 4 3 0 ... $ X4: num 2 3 2 4 1 2 3 3 3 1 ... $ X5: num 4 4 3 3 3 4 4 4 3 3 ... $ X6: num 4 4 4 3 3 4 4 4 3 3 ... $ X7: num 4 2 3 4 3 3 4 4 3 3 ... ### RSM(test) Warning in sqrt(diag(solve(parest$hessian))) : NaNs produced Warning in sqrt(diag(lres$W %*% solve(parest$hessian) %*% t(lres$W))) : NaNs produced Results of RSM estimation: Call: RSM(X = test) Conditional log-likelihood: 13768080 Number of iterations: 5 Number of parameters: 9 Item (Category) Difficulty Parameters (eta): X2X3 X4X5X6X7 Cat 2 Cat 3 Cat 4 Estimate -700.9534 -456.4338 901.4649 -1322.033 -1256.828 -673.2412 -231.5791 -3979.953 1911.748203 Std.ErrNaN NaN NaN NaN NaN NaN NaN NaN1.036215 ### So, I use a shortened data, run like this: RSM(test[1:283,]) Results of RSM estimation: Call: RSM(X = test[1:283, ]) Conditional log-likelihood: -1015.388 Number of iterations: 29 Number of parameters: 9 Item (Category) Difficulty Parameters (eta): X2 X3X4 X5 X6 X7 Cat 2 Cat 3Cat 4 Estimate -0.0739008 0.05655997 1.2203061 -1.1212188 -0.7547957 -0.0378593 1.2020940 3.2657449 8.501760 Std.Err 0.1018267 0.10026033 0.1003035 0.1195827 0.1125677 0.1013739 0.4387487 0.8161754 1.226057 ### However, when I add one record, from 283 to 284, RSM(test[1:284,]) Warning in sqrt(diag(solve(parest$hessian))) : NaNs produced Warning in sqrt(diag(lres$W %*% solve(parest$hessian) %*% t(lres$W))) : NaNs produced Results of RSM estimation: Call: RSM(X = test[1:284, ]) Conditional log-likelihood: 3369344 Number of iterations: 6 Number of parameters: 9 Item (Category) Difficulty Parameters (eta): X2X3 X4X5X6X7 Cat 2 Cat 3 Cat 4 Estimate -670.9137 -554.6215 527.789997 -1359.721 -1127.137 -626.1859 -126.7199 -4753.367 2134.408482 Std.Err 2152.9823 1679.2875 2.377241 NaN NaN 3394.5180 562.6800 11045.9353.526018 I can’t find any special values in the data list X1 X2 X3 X4 X5 X6 X7 270 4 4 4 3 4 4 4 271 3 3 3 1 3 3 3 272 4 3 4 1 4 4 4 273 3 3 3 3 3 3 3 274 3 3 4 4 3 3 4 275 3 3 3 3 3 3 3 276 1 3 3 3 3 3 3 277 3 3 3 2 4 4 4 278 2 3 2 1 3 3 3 279 3 3 3 2 3 3 3 280 3 3 2 2 3 3 3 281 3 4 4 3 4 3 3 282 2 2 2 2 2 2 2 283 3 4 3 1 4 3 3 284 2 4 2 1 4 3 2 285 3 3 3 3 3 3 3 286 4 4 3 4 4 4 4 287 4 3 4 0 3 4 4 288 0 3 4 0 4 4 1 289 4 4 4 4 4 4 4 290 3 3 3 3 3 3 2 If I input many different numbers of data, the results often become strange. The data is appended at the letter. No na values in all data. Great thanks all data is in the appended file __ 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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] What are the pros and cons of the log.p parameter in (p|q)norm and similar?
Sent off-list Thanks Bill and Duncan. I only asked for advice but I got an education too. Michael On 03/08/2021 21:24, Bill Dunlap wrote: In maximum likelihood problems, even when the individual density values are fairly far from zero, their product may underflow to zero. Optimizers have problems when there is a large flat area. > q <- runif(n=1000, -0.1, +0.1) > prod(dnorm(q)) [1] 0 > sum(dnorm(q, log=TRUE)) [1] -920.6556 A more minor advantage for some probability-related functions is speed. E.g., dnorm(log=TRUE,...) does not need to evaluate exp(). > q <- runif(1e6, -10, 10) > system.time(for(i in 1:100)dnorm(q, log=FALSE)) user system elapsed 9.13 0.11 9.23 > system.time(for(i in 1:100)dnorm(q, log=TRUE)) user system elapsed 4.60 0.19 4.78 -Bill On Tue, Aug 3, 2021 at 11:53 AM Duncan Murdoch <mailto:murdoch.dun...@gmail.com>> wrote: On 03/08/2021 12:20 p.m., Michael Dewey wrote: > Short version > > Apart from the ability to work with values of p too small to be of much > practical use what are the advantages and disadvantages of setting this > to TRUE? > > Longer version > > I am contemplating upgrading various functions in one of my packages to > use this and as far as I can see it would only have the advantage of > allowing people to use very small p-values but before I go ahead have I > missed anything? I am most concerned with negatives but if there is any > other advantage I would mention that in the vignette. I am not concerned > about speed or the extra effort in coding and expanding the documentation. > These are often needed in likelihood problems. In just about any problem where the normal density shows up in the likelihood, you're better off working with the log likelihood and setting log = TRUE in dnorm, because sometimes you want to evaluate the likelihood very far from its mode. The same sort of thing happens with pnorm for similar reasons. Some likelihoods involve normal integrals and will need it. I can't think of an example for qnorm off the top of my head, but I imagine there are some: maybe involving simulation way out in the tails. The main negative about using logs is that they aren't always needed. Duncan Murdoch __ R-help@r-project.org <mailto:R-help@r-project.org> mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help <https://stat.ethz.ch/mailman/listinfo/r-help> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html <http://www.R-project.org/posting-guide.html> and provide commented, minimal, self-contained, reproducible code. <http://www.avg.com/email-signature?utm_medium=email_source=link_campaign=sig-email_content=emailclient> Virus-free. www.avg.com <http://www.avg.com/email-signature?utm_medium=email_source=link_campaign=sig-email_content=emailclient> <#DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2> -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] What are the pros and cons of the log.p parameter in (p|q)norm and similar?
Short version Apart from the ability to work with values of p too small to be of much practical use what are the advantages and disadvantages of setting this to TRUE? Longer version I am contemplating upgrading various functions in one of my packages to use this and as far as I can see it would only have the advantage of allowing people to use very small p-values but before I go ahead have I missed anything? I am most concerned with negatives but if there is any other advantage I would mention that in the vignette. I am not concerned about speed or the extra effort in coding and expanding the documentation. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] if statement and for loop question
Dear Kai
When you ask again it is best to tell us what your input is and what
output you were hoping for and what you actually got. If you can make a
small data-set which shows all that then your post will be much more
likely to get a helpful response. If you want to transfer the data-set
to us then using dput() will ensure that we have exactly the same as you
have. Do not forget to tell us what libraries, if any, you have loaded too.
Michael
On 31/05/2021 17:26, Kai Yang via R-help wrote:
Hi Jim,
Sorry to post "same" question, because
1. I was asking to use plain text format. I have to post my question again. But
I don't know if it is working.
2. I'm a beginner for R (< 2 month). It may not easy for me to ask a "clear" R
question. My current work is to transfer my SAS code into R, especially for data
manipulation part.
I'll do my best to ask a non."same" question later.
Thanks,
Kai
On Sunday, May 30, 2021, 10:44:41 PM PDT, Jim Lemon
wrote:
Hi Kai,
You seem to be asking the same question again and again. This does not
give us the warm feeling that you know what you want.
testdf<-data.frame(a=c("Negative","Positive","Neutral","Random","VUS"),
b=c("No","Yes","No","Maybe","Yes"),
c=c("Off","On","Off","Off","On"),
d=c("Bad","Good","Bad","Bad","Good"),
stringsAsFactors=FALSE)
testdf
match_strings<-c("Positive","VUS")
testdf$b<-ifelse(testdf$a %in% match_strings,testdf$b,"")
testdf$c<-ifelse(testdf$a %in% match_strings,testdf$c,"")
testdf$d<-ifelse(testdf$a %in% match_strings,testdf$d,"")
testdf
I have assumed that you mean "zero length strings" rather than
"zeros". Also note that your initial code was producing logical values
that were never assigned to anything.
Jim
On Mon, May 31, 2021 at 2:29 AM Kai Yang via R-help
wrote:
Hello List,I have a data frame which having the character columns:
| a1 | b1 | c1 | d1 |
| a2 | b2 | c2 | d2 |
| a3 | b3 | c3 | d3 |
| a4 | b4 | c4 | d4 |
| a5 | b5 | c5 | d5 |
I need to do: if a1 not = "Positive" and not = "VUS" then values of b1, c1 and
d1 will be zero out. And do the same thing for the a2 to a5 series.
I write the code below to do this. But it doesn't work. Would you please
correct my code?
Thank you,
Kai
for (i in 1:5)
{
if (isTRUE(try$a[i] != "Positive" && try$a[i] != "VUS"))
{
try$b[i]== ''
try$c[i] == ''
try$d[i]== ''
}
}
[[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.
[[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.
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] calculating area of ellipse
Dear Stephen
In that application the axes would be sensitivity and specificity (or
their inverses) or some transformation of them like logits so the units
would be the same. Whether the area has any scientific meaning I am not
sure.
Michael
On 11/05/2021 15:20, Stephen Ellison wrote:
In doing meta-analysis of diagnostic accuracy I produce ellipses of confidence
and prediction intervals in two dimensions. How can I calculate the area of
the ellipse in ggplot2 or base R?
There are established formulae for ellipse area, but I am curious: in a 2-d
ellipse with different quantities (eg coefficients for salary and age)
represented by the different dimensions, what does 'area' mean?
S
***
This email and any attachments are confidential. Any u...{{dropped:13}}
__
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] mgcv: bam warning messages and non-convergence
I have a large dataset of 118225 observations from 16 columns and as such I’ve been using bam, rather than gam, for my analyses. The response variable I’m using is count data but it’s overdispersed, and as such, I thought I’d use a negative binomial model. I have 5 explanatory variables, which are biologically important. Two are numerical and 3 are categorical. I’ve only applied a smoother to the first numerical explanatory variable, because, from some prior analyses I found that TL had edf values of 1.01 and was therefore linear. I also have included categorical two random effects in the model. m3 <- bam(deg ~ s(SE_score) + TL + species + sex + season + year + s(code, bs = 're') + s(monthyear, bs = 're'), family=nb(), data=node_dat, method = "REML") th <- m3$family$getTheta(TRUE) #extracts theta m3 <- bam(deg ~ s(SE_score) + TL + species + sex + season + year + s(code, bs = 're') + s(monthyear, bs = 're'), family=nb(th), data=node_dat, method = "REML") summary(m3) However I’m getting this warning and I can’t find out what it means There were 32 warnings (use warnings() to see them) > warnings() Warning messages: 1: In pmax(1, y)/mu : longer object length is not a multiple of shorter object length 2: In y * log(pmax(1, y)/mu) : longer object length is not a multiple of shorter object length Is this an issue? The model converges, and I’ve checked overdispersion again and get this value > E3 <- resid(m3, type = "pearson") > sum(E3^2)/m3$df.res [1] 0.7436045 So this suggests there is some under dispersion now? Also the model summary gives > summary(m3) Family: Negative Binomial(0.055) Link function: log I’ve read that the 0.055 is also a measure of dispersion so which one is correct? I was confused about all this and I have a lot of zeros in my data (about 96%) so I thought I’d also try an zero inflated poisson, however is does not converge. m4 <- bam(deg ~ s(SE_score) + TL + species + sex + season + year + s(code, bs = 're') + s(monthyear, bs = 're'), family=ziP(), data=node_dat, method = "REML") Warning message: In bgam.fit(G, mf, chunk.size, gp, scale, gamma, method = method, : algorithm did not converge Is there any reason why it does not onverge? And maybe a zero inflated negative binomial would better but I’m not sure how to undertake that. I know there’s a lot here but any help would be appreciated. Many thanks, Mike Michael Williamson London NERC DTP Candidate Email: michael.william...@kcl.ac.uk<mailto:michael.william...@kcl.ac.uk> Phone: +447764836592 Skype: mikejwilliamson Twitter: @mjw_marine Website: www.thenetlab.uk<https://eur03.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.thenetlab.uk%2F=01%7C01%7Cmichael.williamson%40kcl.ac.uk%7C07c592826b364b249c9208d84e5dbc12%7C8370cf1416f34c16b83c724071654356%7C0=vaibGznfTGGiS7l0lHuRaQ3w4fnEQGaXIfgQ34OrhG4%3D=0> Most recent paper: Williamson, M. J. et al. (2021). Analysing detection gaps in acoustic telemetry data to infer differential movement patterns in fish. Ecology and Evolution, 11, 2717-2730. https://doi.org/10.1002/ece3.7226 [[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] readline in function call with space in prompt.
The function test as defined below by Jeremie works as I would have
expected for me on Windows so I am unable to replicate the problem there.
R version 4.0.3 (2020-10-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19041)
Matrix products: default
locale:
[1] LC_COLLATE=English_United Kingdom.1252
[2] LC_CTYPE=English_United Kingdom.1252
[3] LC_MONETARY=English_United Kingdom.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United Kingdom.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] compiler_4.0.3
On 09/02/2021 09:37, Martin Maechler wrote:
Jeremie Juste
on Mon, 08 Feb 2021 14:28:33 +0100 writes:
> Hello,
> I have noticed a behavior that I don't understand. When I call the
> following function from the prompt.
> test <- function(){
> a <- readline("selection: ")
> a
> }
>> test()
>> selection: |
> I can only type one character and the readline function exits before I
can
> press enter.
> however
> test1 <- function(){
> a <- readline("selection:")
> a
> }
>> test1()
>> selection:|
> works as expected.
>> selection: abc[Ret]
> However calling directly readline with a space in the prompt does what I
> would expect.
>> a <- readline("selection: ")
>> selection: abc[Ret]
>> a
>> "abc"
> It is the expected behavior or am I missing something?
> Best regards,
> Jeremie
> --
> Jeremie Juste
>> R version 4.0.3 (2020-10-10)
Given that the above works fine in Linux (for Jim Lemon and Rolf Turner),
could you tell us *how* you use R?
In the (Windows) RGui or from Rstudio or ESS or yet another way?
Usually the UI (user interface) should not matter, but rather
the R version etc.
But the UI may be important for a function like readline()
which does UI ..
Martin
__
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.
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] make new collumns with conditions
Dear Krissie
I think you misunderstood Rui's response. He was generating some fake
data to test the code not suggesting you rebuild your data frame.
Michael
On 25/01/2021 16:01, krissievdh wrote:
Hi,
Thanks for your response.
I do get what you're doing. However, the table I sent is just a small piece
of the complete database. So for me to have to add in everything with
structure list (c ..) by hand would be too much work.
Just to give you an idea, the database is around 16000 rows and has 40
columns with other variables that I do want to keep. So I kind of want to
find a way to keep everything and just add a couple of columns with the
calculated time for vigilant behavior and the percentage.
Still thanks for thinking with me. I am looking into the aggregate
function. Hopefully, this could be a solution.
krissie
Op ma 25 jan. 2021 16:44 schreef Rui Barradas :
Hello,
Try the following.
First aggregate the data, then get the totals, then the percentages.
Finally, put the species in the result.
agg <- aggregate(formula = `duration(s)` ~ `observation nr` + `behavior
type`,
data = d_vigi,
FUN = sum,
subset = `behavior type` == 'Vigilant')
agg$total <- tapply(d_vigi$`duration(s)`, d_vigi$`observation nr`, FUN =
sum)
agg$percent <- round(100 * agg$`duration(s)`/agg$total)
res <- merge(agg, d_vigi[c(1, 3:4)])
res[!duplicated(res), ]
Data in dput format:
d_vigi <-
structure(list(`behavior type` = c("Non-vigilant", "Vigilant",
"Vigilant", "Non-vigilant", "Vigilant", "Vigilant", "Non-vigilant",
"Unkown"), `duration(s)` = c(5L, 2L, 2L, 3L, 7L, 2L, 1L, 2L),
`observation nr` = c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L), species =
c("red deer",
"red deer", "red deer", "red deer", "red deer", "red deer",
"red deer", "red deer")), class = "data.frame", row.names = c(NA,
-8L))
Hope this helps,
Rui Barradas
Às 13:57 de 25/01/21, krissievdh escreveu:
Hi,
I have a dataset (d_vigi)with this kind of data:
behavior type duration(s) observation nr species
Non-vigilant 5 1 red deer
Vigilant 2 1 red deer
Vigilant 2 1 red deer
Non-vigilant 3 1 red deer
Vigilant 7 2 red deer
Vigilant 2 2 red deer
Non-vigilant 1 2 red deer
Unkown 2 2 red deer
Now I have to calculate the percentage of vigilant behavior spent per
observation.
So eventually I will need to end up with something like this:
Observation nr Species vigilant(s) total (s) percentage of vigilant (%)
1 red deer 4 12 33
2 red deer 9 12 75
Now I know how to calculate the total amount of seconds per observation.
But I don't know how I get to the total seconds of vigilant behavior per
observation (red numbers). If I could get there I will know how to
calculate the percentage.
I calculated the total duration per observation this way:
for(id in d_vigi$Obs.nr){
d_vigi$t.duration[d_vigi$Obs.nr==id]<-sum(d_vigi$'Duration.(s).x'[d_vigi$Obs.nr==id])
}
this does work and gives me the total (s) but i don't know how to get to
the sum of the seconds just for the vigilant per observation number. Is
there anyone who could help me?
Thanks,
Krissie
[[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.
[[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.
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] make new collumns with conditions
Dear Krissie
I think you may be looking for the aggregate command.
Note that this is a plain text list so if you post in HTML we do not see
what you see. In this case we did not see any red numbers.
Michael
On 25/01/2021 13:57, krissievdh wrote:
Hi,
I have a dataset (d_vigi)with this kind of data:
behavior type duration(s) observation nr species
Non-vigilant 5 1 red deer
Vigilant 2 1 red deer
Vigilant 2 1 red deer
Non-vigilant 3 1 red deer
Vigilant 7 2 red deer
Vigilant 2 2 red deer
Non-vigilant 1 2 red deer
Unkown 2 2 red deer
Now I have to calculate the percentage of vigilant behavior spent per
observation.
So eventually I will need to end up with something like this:
Observation nr Species vigilant(s) total (s) percentage of vigilant (%)
1 red deer 4 12 33
2 red deer 9 12 75
Now I know how to calculate the total amount of seconds per observation.
But I don't know how I get to the total seconds of vigilant behavior per
observation (red numbers). If I could get there I will know how to
calculate the percentage.
I calculated the total duration per observation this way:
for(id in d_vigi$Obs.nr){
d_vigi$t.duration[d_vigi$Obs.nr==id]<-sum(d_vigi$'Duration.(s).x'[d_vigi$Obs.nr==id])
}
this does work and gives me the total (s) but i don't know how to get to
the sum of the seconds just for the vigilant per observation number. Is
there anyone who could help me?
Thanks,
Krissie
[[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.
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] [SPAM] Re: Different results on running Wilcoxon Rank Sum test in R and SPSS
See comments inline
On 19/01/2021 10:46, bharat rawlley wrote:
Thank you for the reply and suggestion, Michael!
I used dput() and this is the output I can share with you. Simply
explained, I have 3 columns namely, drug_code, freq4w_n and PFD_n. Each
column has 132 values (including NA). The problem with the Wilcoxon Rank
Sum test has been described in my first email.
Please do let me know if you need any further clarification from my
side! Thanks a lot for your time!
structure(list(drug_code = c(0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0,
1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1,
0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1,
1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1,
0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0,
1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1,
1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 0), freq4w_n = c(1,
NA, NA, 0, NA, 4, NA, 10, NA, 0, 6, NA, NA, NA, NA, NA, 10, NA,
0, NA, NA, NA, NA, 0, NA, 0, NA, NA, NA, 0, NA, 0, NA, NA, NA,
NA, NA, NA, NA, NA, 0, 0, 12, 0, NA, 1, 2, 1, 2, 2, NA, 28, 0,
NA, 4, NA, 1, NA, NA, NA, NA, NA, 0, 3, 1, NA, NA, NA, NA, 4,
28, NA, NA, 0, 2, 12, 0, NA, NA, NA, 0, NA, 0, NA, NA, NA, NA,
NA, NA, NA, NA, NA, 3, NA, NA, NA, NA, NA, NA, 6, 1, NA, NA,
NA, 0, NA, NA, NA, 0, 0, NA, 0, NA, 2, 8, 3, NA, NA, NA, 0, NA,
NA, NA, 9, NA, NA, NA, NA, NA, NA, NA, NA), PFD_n = c(27, NA,
NA, 28, NA, 26, NA, 20, NA, 30, 24, NA, NA, NA, NA, NA, 18, NA,
28, NA, NA, NA, NA, 28, NA, 28, NA, NA, NA, 28, NA, 28, NA, NA,
NA, NA, NA, NA, NA, NA, 28, 28, 16, 28, NA, 27, 26, 27, 26, 26,
NA, 0, 30, NA, 24, NA, 27, NA, NA, NA, NA, NA, 28, 25, 27, NA,
NA, NA, NA, 26, 0, NA, NA, 28, 26, 16, 28, NA, NA, NA, 28, NA,
28, NA, NA, NA, NA, NA, NA, NA, NA, NA, 25, NA, NA, NA, NA, NA,
NA, 22, 27, NA, NA, NA, 28, NA, NA, NA, 28, 28, NA, 28, NA, 26,
20, 25, NA, NA, NA, 30, NA, NA, NA, 19, NA, NA, NA, NA, NA, NA,
NA, NA)), row.names = c(NA, -132L), class = c("tbl_df", "tbl",
"data.frame"))
Yours sincerely
Bharat Rawlley
On Tuesday, 19 January, 2021, 03:53:27 pm IST, Michael Dewey
wrote:
Unfortunately your data did not come through. Try using dput() and then
pasting that into the body of your e-mail message.
On 18/01/2021 17:26, bharat rawlley via R-help wrote:
> Hello,
> On running the Wilcoxon Rank Sum test in R and SPSS, I am getting the
following discrepancies which I am unable to explain.
> Q1 In the attached data set, I was trying to compare freq4w_n in
those with drug_code 0 vs 1. SPSS gives a P value 0.031 vs R gives a P
value 0.001779.
> The code I used in R is as follows -
> wilcox.test(freq4w_n, drug_code, conf.int = T)
If I store your data in dat and then go
wilcox.test(freq4w_n ~ drug_code, dat)
I get a p-value of 0.031 agreeing with SPSS
The reason you are getting something different is that you are not
specifying the first two parameters to wilcox.test() correctly.
>
>
> Q2 Similarly, in the same data set, when trying to compare PFD_n in
those with drug_code 0 vs 1, SPSS gives a P value 0.038 vs R gives a P
value < 2.2e-16.
> The code I used in R is as follows -
> wilcox.test(PFD_n, drug_code, mu = 0, alternative = "two.sided",
correct = TRUE, paired = FALSE, conf.int = TRUE)
>
>
> I have tried searching on Google and watching some Youtube tutorials,
I cannot find an answer, Any help will be really appreciated, Thank you!
> __
> R-help@r-project.org <mailto:R-help@r-project.org> mailing list -- To
UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
<https://stat.ethz.ch/mailman/listinfo/r-help>
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
<http://www.R-project.org/posting-guide.html>
> and provide commented, minimal, self-contained, reproducible code.
>
--
Michael
http://www.dewey.myzen.co.uk/home.html
<http://www.dewey.myzen.co.uk/home.html>
<http://www.avg.com/email-signature?utm_medium=email_source=link_campaign=sig-email_content=emailclient>
Virus-free. www.avg.com
<http://www.avg.com/email-signature?utm_medium=email_source=link_campaign=sig-email_content=emailclient>
<#DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2>
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] Different results on running Wilcoxon Rank Sum test in R and SPSS
Unfortunately your data did not come through. Try using dput() and then pasting that into the body of your e-mail message. On 18/01/2021 17:26, bharat rawlley via R-help wrote: Hello, On running the Wilcoxon Rank Sum test in R and SPSS, I am getting the following discrepancies which I am unable to explain. Q1 In the attached data set, I was trying to compare freq4w_n in those with drug_code 0 vs 1. SPSS gives a P value 0.031 vs R gives a P value 0.001779. The code I used in R is as follows - wilcox.test(freq4w_n, drug_code, conf.int = T) Q2 Similarly, in the same data set, when trying to compare PFD_n in those with drug_code 0 vs 1, SPSS gives a P value 0.038 vs R gives a P value < 2.2e-16. The code I used in R is as follows - wilcox.test(PFD_n, drug_code, mu = 0, alternative = "two.sided", correct = TRUE, paired = FALSE, conf.int = TRUE) I have tried searching on Google and watching some Youtube tutorials, I cannot find an answer, Any help will be really appreciated, Thank you! __ 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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] Help with connection issue for R (just joined, leading R for our agency)
Just to add to Petr's comment there are other basic editors with syntax highlighting like Notepad++ which are also OK if you want a fairly minimalist approach. Michael On 14/12/2020 08:16, PIKAL Petr wrote: Hallo Alejandra Although RStudio and ESS could help with some automation (each with its own way), using R alone is not a big problem, especially if you are not familiar with Emacs basics and perceiving RStudio issues. I use R with simple external editor - it could be notepad but I could recommend TINN-R https://sourceforge.net/projects/tinn-r/ which has syntax highlighting and works smoothly if R is console is set to multiple windows. It is also quite easy to manage. Good luck with R. Cheers Petr -Original Message- From: R-help On Behalf Of Alejandra Barrio Gorski Sent: Tuesday, December 8, 2020 7:48 PM To: R-help@r-project.org Subject: [R] Help with connection issue for R (just joined, leading R for our agency) Dear fellow R users, Greetings, I am new to this list. I joined because I am pioneering the use of R for the agency I work for. I essentially work alone and would like to reach out for help on an issue I have been having. Here it is: - From one day to the next, my RStudio does not execute commands when I press ctrl + enter. Nothing happens, and then after a few minutes out of nowhere, it runs everything at once. This makes it very hard to do my work. - I tried uninstalling and re-installing both R and Rstudio, but the error comes up again. I tested commands on my R program alone, and it works fine there. It could be the way that Rstudio connects to R. - I am on a Windows 10 computer. I work for a government agency so there may be a few firewall/virus protection issues. I would love any pointers. Thank you, Alejandra -- *Alejandra Barrio* Linkedin <https://www.linkedin.com/in/alejandra-barrio/> | Website <https://www.ocf.berkeley.edu/~alejandrabarrio/> MPP | M.A., International and Area Studies University of California, Berkeley [[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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] Find the ideal cluster
Dear Jovani
If you cross-post on CrossValidated as well as here it is polite to give
a link so people do not answer here when someone has already answered
there, or vice versa.
Michael
On 12/12/2020 15:27, Jovani T. de Souza wrote:
So, I and some other colleagues developed a hierarchical clustering
algorithm to basically find the main clusters involving agricultural
industries according to a particular city (e.g. London city).. We
structured this algorithm in R. It is working perfectly. So, according to
our filters that we inserted in the algorithm, we were able to generate 6
clustering scenarios to London city. For example, the first scenario
generated 2 clusters, the second scenario 5 clusters, and so on. I would
therefore like some help on how I can choose the most appropriate one. I
saw that there are some packages that help in this process, like `pvclust`,
but I couldn't use it for my case. I am inserting a brief executable code
below to show the essence of what I want.
Any help is welcome! If you know how to use using another package, feel
free to describe.
Best Regards.
library(rdist)
library(geosphere)
library(fpc)
df<-structure(list(Industries = c(1,2,3,4,5,6),
+Latitude = c(-23.8, -23.8, -23.9, -23.7,
-23.7,-23.7),
+Longitude = c(-49.5, -49.6, -49.7, -49.8,
-49.6,-49.9),
+Waste = c(526, 350, 526, 469, 534, 346)), class =
"data.frame", row.names = c(NA, -6L))
df1<-df
#clusters
coordinates<-df[c("Latitude","Longitude")]
d<-as.dist(distm(coordinates[,2:1]))
fit.average<-hclust(d,method="average")
clusters<-cutree(fit.average, k=2)
df$cluster <- clusters
> df
Industries Latitude Longitude Waste cluster
1 1-23.8 -49.5 526 1
2 2-23.8 -49.6 350 1
3 3-23.9 -49.7 526 1
4 4-23.7 -49.8 469 2
5 5-23.7 -49.6 534 1
6 6-23.7 -49.9 346 2
>
clusters1<-cutree(fit.average, k=5)
df1$cluster <- clusters1
> df1
Industries Latitude Longitude Waste cluster
1 1-23.8 -49.5 526 1
2 2-23.8 -49.6 350 1
3 3-23.9 -49.7 526 2
4 4-23.7 -49.8 469 3
5 5-23.7 -49.6 534 4
6 6-23.7 -49.9 346 5
>
[[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.
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] scott-knot ESD effect size test
There seems to be a package on CRAN dedicated exclusively to this test. It is called ScottKnotESD rather unoriginally. Michael. On 02/12/2020 10:32, Neha gupta wrote: I have the following data from resample svm= svm$resample$RMSE nn= nn$resample$RMSE we perform the statistical tests like wilcox.test(svm, nn) I have a question, can we perform the scott-knot ESD test here? if yes, how? [[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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] WG: Packegs for mathematics
Dear Mihaela As well as the search tips you have been given you may want to look at the CRAN Task Views. There is one for Multivariate https://cran.r-project.org/view=Multivariate and one for NumericalMathematics https://cran.r-project.org/view=NumericalMathematics and several others which may be relevant. Note that the task views tell you what exists but they do not offer best buy information. For that you may need to come back here and ask about a specific subset of packages. Michael On 16/11/2020 10:02, ELLENBERGER Mihaela via R-help wrote: Von: ELLENBERGER Mihaela Gesendet: Montag, 16. November 2020 10:24 An: cran-submissi...@r-project.org Betreff: Packegs for mathematics Hello I'm undergraduated student of Biomedical Sciences in Switzerland and learning how to work with R. Currently I'm using R-Studio. I would like to now, which packeges sholud I install please for mathematics ( analysis)? My collegeues don't use all the same maths software ( some of them are using Excel, SPSS, Mathematica, etc.) -Funktionen,Serien, Folgen,Grenzwerte, Stetigkeiten -Integralrechnungen -Differentialrechnungen -Analysis multivariant Could you sent me an advice please, because there are so many packeges in R. Kind regards Mihaela Ellenberger [[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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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 using "revtools" package
Dear John
Your .bib file did not make it through the system as only a very limited
set of attachments is supported. Try to make a minimal .bib file which
triggers the problem by splitting your file in half, testing each half,
repeat until you get the smallest file which triggers the error.
Hopefully that will be a single entry which you can then include in your
post. It also increases the chances of someone looking at it as a large
.bib file is going to deter people.
Michael
On 05/11/2020 02:38, John wrote:
Hi,
I tried to read a bib file by "read_bibliography" function in
"revtools" package, but I got an error message but don't know how to fix
it. Anyone can help? Thank you so much!!
###
library(revtools)
data2 <- read_bibliography("mlf_ref201105_2.bib", return_df = FALSE)
###
Error message:
Error in if (any(names(entry_list) == "author")) { :
missing value where TRUE/FALSE needed
My "mlf_ref201105_2.bib" file is attached.
I have no problem reading the built-in ris file:
file_location <- system.file(
"extdata",
"avian_ecology_bibliography.ris",
package = "revtools")
x <- read_bibliography(file_location)
__
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.
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] build a literature database
Dear John If you are doing a systematic review have you thought of investigating some of the packages in the CRAN Task View on MetaAnalysis like metagear or revtools? I have not used any of them but they claim to support that part of the process. Michael On 04/11/2020 09:22, John wrote: Hi, I 'd like to create a table for literature review. Is there any good data structure (database) I may use? Now I just use a simple dataframe as follows, but in the second item I swap the order of year and author, and it records 2013 as author and "XH" as year. Is there any better structure ? If not, how may I fix the current condition?Thanks! df <- data.frame(author = NA, year = NA, title = NA, country = NA, sample = NA, data = NA, result = NA, note = NA, stringsAsFactors=FALSE) df[1, ]<-c( author = "Moore", year = 2020, title = "Statistics and data analysis", country = "Colombia", sample = NA, data = "firm level", result = NA, note = NA) df[nrow(df)+1,]<- c(year = 2013, author = "XH", title = NA, country = NA, sample = NA, data = NA, result = NA, note = NA) [[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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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: variable not found
Dear Hannah I think the problem is that attach() is not doing what you think it is. It does seem to make it easy to make mistakes. I would suggest switching to using with() instead or using the data = parameter to functions which support it. Michael On 30/10/2020 08:15, Hannah Van Impe wrote: Hello I have a question. I made an r-script and did a few commands needed to make some new variables. They all work out well, and when I run the commands, the new variables appear in the dataset. I can also work with these new variables to make other new variables from them. Also, when I use summary(dataset), those new variables appear in the summary of the dataset. But, when I do summary(new variable) => error: variable not found. For example: summary(couple_id) => Error in summary(couple_id) : object 'couple_id' not found. What can I do about this? R-script: attach(ipumsi_8_dta) library(tinytex) library(dplyr) library(ggplot2) library(tidyr) library(knitr) library(forcats) library(mice) library(pander) library(ggcorrplot) library(lubridate) # true/false code when sploc is greater than zero and sprule is equal to 1 or 2 ipumsi_8_dta <- mutate(ipumsi_8_dta, rule_union = sploc>0 & (sprule==1 | sprule==2)) ## creating numeric code for rule_union & rule_unionn: 1 when sploc is greater than zero and sprule is equal to 1 or 2, 0 if not. ## This is neccesary because otherwise it is a logical code and we cannot multiply with it, which is needed ipumsi_8_dta <- ipumsi_8_dta %>% mutate(rule_union_numeric = as.numeric(rule_union)) ### creating unique numeric code for sploc / pernum variables ipumsi_8_dta <- mutate(ipumsi_8_dta, sploc_pernum_code = pernum*(sex==1) + sploc*(sex==2)) dividing serial by 1000, otherwise, the ultimate couple_id is too large, and it works in this dataset because the serials start at 1000 (I will have to check if this works for other datasets) ipumsi_8_dta <- mutate(ipumsi_8_dta, serial_divided = serial%/%1000) # creating unique union identifier ipumsi_8_dta$union_id = paste0(ipumsi_8_dta$country, ipumsi_8_dta$year, ipumsi_8_dta$serial_divided, ipumsi_8_dta$sploc_pernum_code) ipumsi_8_dta <- ipumsi_8_dta %>% mutate(union_id_numeric = as.numeric(union_id)) ipumsi_8_dta <- mutate(ipumsi_8_dta, couple_id = union_id_numeric*rule_union_numeric) ## Understanding the couple_id variable summary(couple_id) summary(ipumsi_8_dta) I also attached my dataset. Thank you very much for the answer! Hannah Van Impe __ 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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] Need help in programming R on functional data
You have already asked this and people gave you a variety of answers. Just asking again without clarifying why those answers did not help you is not going to solve your problem. Tell us what you tried and why it failed might help. Michael On 29/10/2020 07:50, Ablaye Ngalaba wrote: Hello, I need to generate the functional data in R programming but I can't do it. Please, an example of R code that generates the functional data can help me. Thank you and have a nice day [[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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] Dear R experts, I have a question about meta-analysis
Dear Zhang There is a mailing list dedicated to meta-analysis in R where your question may get more attention. Before posting it would be a good idea to look at the archives as this issue comes up there repeatedly. https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis// For a link to the archive and for instructions on registering for the list. Michael On 27/10/2020 13:32, 21803...@zju.edu.cn wrote: Dear R experts, Greetings from China! I'm Zhang in the College of Education, Zhejiang University, and I am recently running a meta-analysis. Since research using the randomized controlled trial (RCT) often dismissed reporting the correlation (r) between multivariate outcomes, for instance, a study measuring students' gains on problem solving skills with three aspects and reported the pre-post scores respectively, but obviously these three aspects were correlated. I wonder if and how I could integrate the effect sizes among these three aspects into an overall effect size without getting the concrete 'r'? Could R project and (or) function help me solve this problem? Best, Zhang -- Zhang Enming (张恩铭) PhD Student Department of Curriculum and Instruction College of Education, Zhejiang University Tel: +86 17649850218 Email: 21803...@zju.edu.cn [[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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] vanilla session in R Gui or RStudio
[env: Windows, R 3.6.6] When I start R from the R Gui icon or from RStudio, I get a large number of packages loaded via a namespace. Not entirely clear where these come from. As a result, I often run into problems updating packages because something is already loaded. How can start a new gui session with minimal packages loaded? > sessionInfo() R version 3.6.3 (2020-02-29) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 7 x64 (build 7601) Service Pack 1 Matrix products: default locale: [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C [5] LC_TIME=English_United States.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base loaded via a namespace (and not attached): [1] statmod_1.4.34 xfun_0.18tidyselect_1.1.0 reshape2_1.4.4 purrr_0.3.4 mitools_2.4 [7] splines_3.6.3lattice_0.20-41 coefplot_1.2.6 carData_3.0-4 colorspace_1.4-1 vctrs_0.3.4 [13] generics_0.0.2 htmltools_0.5.0 yaml_2.2.1 survival_3.2-7 rlang_0.4.7 pillar_1.4.6 [19] nloptr_1.2.2.2 glue_1.4.2 DBI_1.1.0lifecycle_0.2.0 plyr_1.8.6 stringr_1.4.0 [25] effects_4.2-0munsell_0.5.0gtable_0.3.0 evaluate_0.14 knitr_1.30 fansi_0.4.1 [31] Rcpp_1.0.5 scales_1.1.1 useful_1.2.6 fs_1.4.2 lme4_1.1-23 packrat_0.5.0 [37] ggplot2_3.3.2digest_0.6.25stringi_1.4.6insight_0.9.6 dplyr_1.0.2 survey_4.0 [43] grid_3.6.3 cli_2.1.0tools_3.6.3 magrittr_1.5 tibble_3.0.4 crayon_1.3.4 [49] pkgconfig_2.0.3 ellipsis_0.3.1 MASS_7.3-53 Matrix_1.2-18 reprex_0.3.0 assertthat_0.2.1 [55] minqa_1.2.4 rmarkdown_2.4rstudioapi_0.11 R6_2.4.1 boot_1.3-25 nnet_7.3-14 [61] nlme_3.1-149 compiler_3.6.3 > Michael Friendly Email: friendly AT yorku DOT ca Professor, Psychology Dept. & Former Chair, ASA Statistical Graphics Section York University Voice: 416 736-2100 x66249 4700 Keele StreetWeb: http://www.datavis.ca | @datavisFriendly Toronto, ONT M3J 1P3 CANADA __ 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] R for mac
Hi I am club mentor for a group of high school students learning R. The Vice President of Information Technology is hoping to broaden and deepen her skill set so she can help others. She is doing well in helping students install R on the Windows operating system. However, a few have problems installing R on mac and one student is struggling to install R on Chromebook. Is there someone who would be willing to provide some pointers? Thank you for considering this request, Michael On Thu, Sep 24, 2020 at 9:10 AM Dirk Eddelbuettel wrote: > > On 23 September 2020 at 15:32, Michael Johnston wrote: > | Hi > | I am club mentor for a group of high school students learning R. A few > have > | problems installing R on mac. The Vp of information technology is hoping > | for training so she can help others. Could you recommend someone to help > | her? > | Thank you for considering this request > > I have nothing to do with R on macOS -- but this mailing list is the place > for mac-focussed discussion. Maybe try that, and you likely need to > subscribe before you can post. > > Dirk > > | Michael > > -- > https://dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org > [[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] Stats help for dissertation project
Dear Raija This list is primarily for R programming questions so most of your post is off-topic here. If you are registered for a degree presumably someone is paying a fee to your institution and someone there is being paid to supervise your project so I would have thought they would be the first port of call. If they fail to meet their obligations then there is a site Cross Validated where you may have better luck. Michael On 18/09/2020 18:19, Raija Hallam wrote: Hello, I am a conservation Masters student who is new to R and in need of some confirmation of my methods/ stats help! My dissertation project is looking at the micronutrient intake (iron, zinc, calcium, magnesium and folic acid/folate) of 18 female monkeys, 14 of which are reproductively 'successful' (their infant survived past 1 year) and 4 of which are reproductively 'unsuccessful' (their infant did not survive to 1 year). The females go through a number of stages throughout their pregnancy, and I would like to focus on 2 of these stages, early gestation and early lactation, as these appear in the literature to be important stages in terms of nutrition. Each stage is broken down into 4 week periods so I have G1, G2 and G3 as early gestation and L1, L2, L3 and L4 as early lactation. These could also be combined into just 2 reproductive stages; early gestation (EG) and early lactation (EL), to make the model a bit simpler. *I first would like to investigate how micronutrient intake is affected by the reproductive stage of females.* To investigate this I am thinking of doing a multivariate multiple regression general linear model, controlling for Female ID: ironintake ~ repphase + (1/FemaleID) zincintake ~ repphase + (1/FemaleID) etc. *I would also like to investigate how the micronutrients they intake affect reproductive 'successfulness'.* To investigate this I am thinking of doing a binomial logistic regression generalised linear model, again controlling for Female ID: Repsuccess ~ ironintake + zincintake + calciumintake etc... +(1/FemaleID) As I'm new to R and a bit rusty with my stats knowledge I would be very grateful for any comments on my current stats tests and methods outlined above. If there are other tests that would fit my data better then I'm all ears! Thank you! [[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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] question including crossover trials in meta-analysis
Dear Vera In addition to what you already have you might like to know about the mailing list specifically dedicated to meta-analysis in R. https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis// You might like to search the archives first as this sort of issue does come up there. You are also likely to get more immediate help if you frame your question in terms of either meta or metafor the two packages which dmetar uses. The authors of meta and metafor all frequent the list Michael On 14/09/2020 20:41, Marc Schwartz via R-help wrote: Hi, Bert has pointed you to some R specific packages for meta-analyses via the Task View. It sounds like you may need to first address some underlying conceptual issues, which strictly speaking, is off-topic for this list. That being said, a quick Google search came up with some possible resources, beyond the Cochrane reference: https://academic.oup.com/ije/article/31/1/140/655940 https://onlinelibrary.wiley.com/doi/abs/10.1002/jrsm.1236 https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0133023 Perhaps with additional conceptual background, that might assist in enabling you to make use of the R packages that provide relevant functionality. Another option, if you want to stay with the dmetar package, would be to contact the package maintainer for some guidance relative to how to use the functionality in the package given your specific use case. Regards, Marc Schwartz On Sep 14, 2020, at 3:14 PM, Bert Gunter wrote: Did you first try a web search? -- you should always do this before posting here. "meta-analysis in R" brought up this: https://CRAN.R-project.org/view=MetaAnalysis Have you looked at this task view yet? 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 Mon, Sep 14, 2020 at 11:50 AM Belgers, V. (Vera) < v.belg...@amsterdamumc.nl> wrote: Dear sir/madam, Thank you in advance for taking the time to read my question. I am currently trying to conduct a meta-analysis combining parallel and crossover trials. According to the Cochrane Handbook, I can include crossover trials by using t-paired statistics. So far, I have managed to conduct a meta-analysis and forest plot of the parallel trials using the dmetar package, but I did not succeed in including the crossover trials. I do have the raw data of most of these crossover trials. Does anybody know how to add crossover trials to the meta-analysis? With kind regards, Vera Belgers __ 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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] truth[(truth[, 1]=="G3" & truth[, 2]=="G2") | (truth[, 1]=="G2" & truth[, 2]=="G3"), 3]<-1
I am afraid this is completely unreadable because you posted in HTML ad
this is a plain text list. Best to resend it having set your mailer to
send plain text as HTML gets mangled here.
Michael
On 06/09/2020 10:58, Hesham A. AL-bukhaiti via R-help wrote:
helloout<-read.csv("outbr.csv")truth<-out[,seq(1,2)]for example :
If row1= G1 and row2=G2 , and row 1 = G2 and row 2= G1,make G3=1 # note G1
and G2 are values from 1 to 2000 #if this happend add to thrid column in truth
1 otherwise add 0 as in statment follow
truth<-cbind(as.character(truth[,1]),as.character(truth[,2])
,as.data.frame(rep(0,,dim(out)[1])));#here just G2 and G3, i want make loop to
cam[are all values from G1 to G2000
truth[(truth[,1]=="G3" & truth[,2]=="G2") | (truth[,1]=="G2" &
truth[,2]=="G3"),3]<-1 ###3(Simply they regulate the other. If element A is in the first
group , and it is related to element B in the second group , and element B also in in the first group , and it is related to
element A(the same element in the first group) in the second group , we write 1 and otherwise 0.
this the distination result:
I want this result
G1 G2 G3 D B 1 B D 1 A D 0 B A 1B C 0A B 1
[[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.
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] R 4.0.2 not running on macOS 10.15.6
Greetings, I am a brand-new user to R and want to install and run this on my iMac running macOS 10.15.6. I successfully downloaded the correct package, installed it with no problems (I was not given any options for Tcl/Tk or Texinfo), and started it up as any other normal app. Some introductory commands like license() and help() worked just fine. But when I went to execute demo(), it just beach-balled. I’ve had to force-quit R a number of times. Yes, I read the Mac installation instructions. Thanks in advance. Mike __ 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] hist from a list
Dear Pedro
Some comments in-line
On 30/07/2020 21:16, Pedro páramo wrote:
Hi all,
I attach my code, the think is I want to make a bar plot the last variable
called "bwchist"
so the X axis are "Accion" and the y axis are "reval" values.
I have prove class(bwchist) and says dataframe but its still a list because
it says me
I have prove to unlist, but it doesnt work
hist(bwchist)
Error in hist.default(bwchist) : 'x' must be numeric
So bwchist is not a numeric variable as hist needs. Aboce you said it is
a data frame but data frames are not numeric.
For future reference your example is way too long for anyone to go
through and try to help you. Try next time to reduce it to the absolute
minimum by removing sections while you still get the error. It is also
easier to get help if you can remove unnecessary packages.
It is also unreadable because you are posting in HTML and that makes the
post unreadable as this is a plain text list.
Michael
Or
barplot(bwchist)
Error in barplot.default(bwchist) : 'height' must be a vector or a matrix
library(PerformanceAnalytics)
library(dplyr)
library(tibble)
library(lubridate)
library(PerformanceAnalytics)
library(quantmod)
library(ggplot2)
library(png)
library(grid)
library(RCurl)
library(tidyquant)
library(timetk)
library(data.table)
Acciona<- tq_get("ANA.MC",from = '2019-12-31',get = "stock.prices")
ACS<- tq_get("ACS.MC",from = '2019-12-31',get = "stock.prices")
Aena<- tq_get("AENA.MC",from = '2019-12-31',get = "stock.prices")
Amadeus<- tq_get("AMS.MC",from = '2019-12-31',get = "stock.prices")
ArcelorMittal<- tq_get("MTS.MC",from = '2019-12-31',get = "stock.prices")
BBVA<- tq_get("BBVA.MC",from = '2019-12-31',get = "stock.prices")
Sabadell<- tq_get("SAB.MC",from = '2019-12-31',get = "stock.prices")
Santander<- tq_get("SAN.MC",from = '2019-12-31',get = "stock.prices")
Bankinter<- tq_get("BKT.MC",from = '2019-12-31',get = "stock.prices")
CaixaBank<- tq_get("CABK.MC",from = '2019-12-31',get = "stock.prices")
Cellnex<- tq_get("CLNX.MC",from = '2019-12-31',get = "stock.prices")
Enagas<- tq_get("ENG.MC",from = '2019-12-31',get = "stock.prices")
ENCE<- tq_get("ENC.MC",from = '2019-12-31',get = "stock.prices")
Endesa<- tq_get("ELE.MC",from = '2019-12-31',get = "stock.prices")
Ferrovial<- tq_get("FER.MC",from = '2019-12-31',get = "stock.prices")
Grifols<- tq_get("GRF.MC",from = '2019-12-31',get = "stock.prices")
Iberdrola<- tq_get("IBE.MC",from = '2019-12-31',get = "stock.prices")
Inditex<- tq_get("ITX.MC",from = '2019-12-31',get = "stock.prices")
Colonial<- tq_get("COL.MC",from = '2019-12-31',get = "stock.prices")
IAG<- tq_get("IAG.MC",from = '2019-12-31',get = "stock.prices")
Mapfre<- tq_get("MAP.MC",from = '2019-12-31',get = "stock.prices")
Melia<- tq_get("MEL.MC",from = '2019-12-31',get = "stock.prices")
Merlin<- tq_get("MRL.MC",from = '2019-12-31',get = "stock.prices")
Naturgy<- tq_get("NTGY.MC",from = '2019-12-31',get = "stock.prices")
REE<- tq_get("REE.MC",from = '2019-12-31',get = "stock.prices")
Repsol<- tq_get("REP.MC",from = '2019-12-31',get = "stock.prices")
SGamesa<- tq_get("SGRE.MC",from = '2019-12-31',get = "stock.prices")
Telefonica<- tq_get("TEF.MC",from = '2019-12-31',get = "stock.prices")
Viscofan<- tq_get("VIS.MC",from = '2019-12-31',get = "stock.prices")
Acerinox<- tq_get("ACX.MC",from = '2019-12-31',get = "stock.prices")
Bankia<- tq_get("BKIA.MC",from = '2019-12-31',get = "stock.prices")
CIE<- tq_get("CIE.MC",from = '2019-12-31',get = "stock.prices")
MasMovil<- tq_get("MAS.MC",from = '2019-12-31',get = "stock.prices")
Almirall<- tq_get("ALM.MC",from = '2019-12-31',get = "stock.prices")
Indra<- tq_get("IDR.MC",from ='2019-12-31',get = "stock.prices")
Indra_daily_returns <- Indra %>%
tq_transmute(select = adjusted, # this specifies which column
to select
mutate_fun = periodReturn, # This specifies what to do
with that column
period = "daily", # This argument calculates Daily
returns
col_rename = "idr_returns") # renames the column
Indra_cum_returns <- Indra_daily_returns %>%
mutate(cr = cumprod(1
Re: [R] Fortune nomination .... Re: Looping through a dataframe
Dear Chris Just send it to the maintainer Achim Zeileis to whom I have cc'ed this in case he is not reading the list at the moment. Michael On 25/07/2020 06:21, Chris Evans wrote: I really don't want to put too many Emails to the list but I had the same reaction to that lovely line ... so now my question is: Is there a way to make fortune nominations other than through the list. I see the package has a long list of illustrious authors but I can't find a side channel or open port for nominations despite a bit of searching. TIA and very best to all, Chris - Original Message - From: "John Kane" To: "David Winsemius" Cc: "R. Help Mailing List" Sent: Saturday, 25 July, 2020 02:55:34 Subject: Re: [R] Fortune nomination Re: Looping through a dataframe Yes I think so. On Fri, 24 Jul 2020 at 20:53, David Winsemius wrote: On 7/21/20 2:31 PM, Jim Lemon wrote: I want to get a total for the number of years of data for each company. When I loop through the data After two one liners using `table`: I'm too lazy to provide a difficult way. Jim __ 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. -- John Kane Kingston ON Canada [[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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] Question about citing R packages in an academic paper
Dear Medhi It is good that you are going to cite the packages properly. I do not think it matters too much whether they are in the article itself or in the supplementary material, the important thing is that they are there. You do not have any legal obligation to cite them, as far as I know although I am not a lawyer, but it is a kindness to the authors. Michael On 22/07/2020 06:26, Dr. Mehdi Dadkhah wrote: Hi, I hope you are doing well! I have a question. I and my colleagues wrote a paper by using R language and its packages. We also used some tutorials. We have words count limitation for our paper. In early version of paper, I cited all packages in the reference list of paper (my paper will be published in an academic journal). After that, we had to keep references to tutorials in the reference list and remove reference to R and its packages (I strictly should fit paper to 2000 words). We created an online appendix and listed all packages which we used. Some packages are related to paper. For example, when i used command 'citation()' in R console for package "tidytext", R informed me that i should cite such package as this: Silge, J., & Robinson, D. (2016). tidytext: Text Mining and Analysis Using Tidy Data Principles in R. JOSS, 1(3). https://doi.org/10.21105/joss.00037 Now, all packages and their related papers are not in the reference list of paper, but they will be available in an online appendix in the journal website. Is such practice OK or should I cite all packages in the reference list? I attached online appendix. We will submit our paper to a journal soon. Should I cite papers which are related to packages in the reference list? if I do not cite, i will face with permission problem regard to usage of package? Many thanks! With best regards, Dr. Mehdi Dadkhah __ 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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] Reading help functions
Dear Pedro If you prefer reading in Spanish which from your name may be the case then if you go to CRAN, click on Documentation | Contributed you will find some documentation in Spanish including several aimed at beginners. There are some in Portuguese as well if I have mis-interpreted your name. Of course ultimately you may prefer to read the documents in English because the program help and vignettes are usually in English and familiarity with the English terminology may help. Michael On 22/07/2020 18:20, Pedro páramo wrote: Hi all, I am trying all time to use ?? help function but most of the time it cost me a lot to understand what they are saying, explaining, there is some manual to step by step how to interpret help guides in R. I hope you can understand me because of my english its not the best also. Many thanks [[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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] output from R to simple html
I can't tell what kind of structure you want to output, but one simple thing that I've done along these lines is to use knitr::kable to output a dataframe in either HTML or LaTeX format (for direct inclusion on a web page or to convert to PDF for other uses). -- Mike On Wed, Jul 8, 2020 at 12:02 PM Marc Roos wrote: > > > > I would like to parse some input to an R script and use its result > output (maybe in json) on a web page. I think this shiny framework is a > bit over kill. What is the simplest to implement this? > > __ > 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] how to calculate odd ratios with R?
Dear Luigi
You could try the epitools package which gives a large number of ways of
doing this. I would have thought that using Wald intervals for the log
odds ration was not optimal with small frequencies.
Michael
On 06/07/2020 14:01, Luigi Marongiu wrote:
Hello,
Is it possible to calculate with a single function the odd ratios?
Now I can use this implement:
```
or <- (De/He)/(Dn/Hn) # Disease exposed, Healthy non-exposed
logo <- log(or)
x <- sqrt(((1/De) + (1/He) + (1/Dn) + (1/Hn)))
lower_ci = exp(logo - 1.96*x)
upper_ci = exp(logo + 1.96*x)
cat("OR:", round(or, 3), "(", round(lower_ci, 3), "-", round(upper_ci, 3), ")",
spe = "")
```
for instance,
```
De <-6
Dn <-3
He <-4
Hn <-5
or <- (De/He)/(Dn/Hn)
logo <- log(or)
x <- sqrt(((1/De) + (1/He) + (1/Dn) + (1/Hn)))
lower_ci = exp(logo - 1.96*x)
upper_ci = exp(logo + 1.96*x)
cat("OR:", round(or, 3), "(", round(lower_ci, 3), "-", round(upper_ci, 3), ")",
spe = "")
OR: 2.5 ( 0.37 - 16.889 )
```
Is there a simple function from some package that can also add a
p-value to this test? Or how can I calculate the p-value on my own?
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] how to summarize results from studies?
Dear Frederik There is also a mailing list dedicated to meta-analysis in R https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis// Michael On 01/07/2020 16:40, Marc Schwartz via R-help wrote: Hi, It sounds like you will want to engage in a meta-analysis. There is a CRAN task view here: https://cran.r-project.org/web/views/MetaAnalysis.html that would be relevant in pointing you to tools in R that can support that approach. That being said, the details of specific methodologies and conceptual assistance would be beyond the scope of this list. You should consider consulting a local statistician for assistance with that, if needed. Regards, Marc Schwartz On Jul 1, 2020, at 11:27 AM, Frederik Feys wrote: Hello everyone I have some studies with results from the same outcome scale. I want to merge them into 1 summarised estimated result and its standard deviation. How do I do that in R? Thank you very much for your help! Frederik __ 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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] Using OpenBLAS on Ubuntu
Or maybe R-devel. On Sun, Jun 28, 2020 at 8:51 PM Jeff Newmiller wrote: > > Wouldn't this question make more sense in r-sig-debian? > > On June 28, 2020 7:54:49 PM PDT, Erin Hodgess wrote: > >Hello! > > > >Hope everyone is doing well. > > > >I'm not sure if this is the correct place to post, but I thought I > >would > >start here. > > > >I have R 4.0.2 on Ubuntu 20.04. I would like to use OpenBLAS with it, > >but > >I'm not sure if I can because I did not compile from source. Is that > >possible, please? > > > >Thanks for any help! > > > >Sincerely, > >Erin > > > >Erin Hodgess, PhD > >mailto: erinm.hodg...@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. > > -- > 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-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 in meta-analysis package Metafor-weights =""
The two packages both define a function forest. When you load the second
package R will have told you that it was masking the definition of
forest from the first package. If you had loaded them in the other order
it would have masked the other one.
In fact it masked seven functions in total in this case as the message
told you.
Michael
On 23/06/2020 16:29, K Amoatwi wrote:
Dear Wolfgang,
Yes! The "metafor::forest(result.md)" works
Thank you very much.
Any reason why "forest(result.md)" was not working?
Kobby
On Tue, Jun 23, 2020 at 11:18 AM Viechtbauer, Wolfgang (SP) <
wolfgang.viechtba...@maastrichtuniversity.nl> wrote:
You have loaded the 'meta' package after 'metafor' and then forest() will
try to use the corresponding function from the meta package and not
metafor. With:
metafor::forest(result.md)
it should work.
Best,
Wolfgang
-Original Message-
From: K Amoatwi [mailto:amoatwi...@gmail.com]
Sent: Tuesday, 23 June, 2020 16:37
To: Viechtbauer, Wolfgang (SP)
Cc: r-help@r-project.org
Subject: Re: [R] Error message in meta-analysis package Metafor-weights
=""
Dear Wolfgang,
I have posted the requested information you asked for.
sessionInfo()
R version 4.0.1 (2020-06-06)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18362)
Matrix products: default
locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United
States.1252
[3] LC_MONETARY=English_United States.1252
LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] meta_4.12-0reshape2_1.4.4 metafor_2.4-0 Matrix_1.2-18
loaded via a namespace (and not attached):
[1] Rcpp_1.0.4.6 lattice_0.20-41MASS_7.3-
51.6 grid_4.0.1
[5] plyr_1.8.6 nlme_3.1-
148 magrittr_1.5 stringi_1.4.6
[9] minqa_1.2.4nloptr_1.2.2.1 boot_1.3-
25splines_4.0.1
[13] statmod_1.4.34 lme4_1.1-
23tools_4.0.1stringr_1.4.0
[17] CompQuadForm_1.4.3 compiler_4.0.1
class(result.md)
[1] "rma.uni" "rma"
Thank you
Kobby
On Tue, Jun 23, 2020 at 3:24 AM Viechtbauer, Wolfgang (SP)
wrote:
Dear Kobby,
Please post the output of sessionInfo() and class(result.md).
Best,
Wolfgang
-Original Message-
From: K Amoatwi [mailto:amoatwi...@gmail.com]
Sent: Monday, 22 June, 2020 22:30
To: Viechtbauer, Wolfgang (SP)
Cc: r-help@r-project.org
Subject: Re: [R] Error message in meta-analysis package Metafor-weights
=""
Hi Wolfgang and All,
I am still practising my meta-analysis with the "Metafor" package, I
tried
to run the code for "Forest plot" and got error message as shown below:
forest(result.md)
forest(result.md)
Error in UseMethod("forest") :
no applicable method for 'forest' applied to an object of class
"c('rma.uni', 'rma')"
Thank you in advance for your support
regards
Kobby
On Tue, Jun 16, 2020 at 12:50 PM Viechtbauer, Wolfgang (SP)
wrote:
Dear Amoatwi,
This way of using the escalc() function has been deprecated. It might be
added back once there is actually any benefit from having this
functionality, but for years it just meant that I had to maintain two
different ways of doing the exact same thing without any additional
benefits.
Best,
Wolfgang
--
Wolfgang Viechtbauer, Ph.D., Statistician | Department of Psychiatry and
Neuropsychology | Maastricht University | P.O. Box 616 (VIJV1) | 6200 MD
Maastricht, The Netherlands | +31 (43) 388-4170 | http://www.wvbauer.com
-Original Message-
From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of K
Amoatwi
Sent: Tuesday, 16 June, 2020 4:50
To: r-help@r-project.org
Subject: [R] Error message in meta-analysis package Metafor-weights =""
Dear All,
I am using the example from one of the tutorial about "Metafor" package
and
"escalc" function, to learn how this package can be applied to do
meta-analysi; the code and the data is directly from the tutorials but
"weights=freq" option in the escalc function is given me error message
This is the code below:
library(metafor) # Load package
#DATASET 1: BCG Vaccine Trials
data(dat.bcg) # BCG meta-analytic dataset
##Formula based Specification
##That is, what if I have multiple rows per study, corresponding to
difference treatment groups?
library(reshape2) # Load package for data reshaping
bcg.long <- melt(dat.bcg[, c("trial", "tpos", "tneg", "cpos", "cneg")],
id
= "trial")
bcg.long$pos <- ifelse(bcg.long$var == "tpos" | bcg.long$var == "cpos",
1,
0)
bcg.long$group <- ifelse(bcg.long$var == "tpos" | bcg.long$var ==
"tneg",
1, 0)
##sample of the data, the first 6 rows
head(bcg.long)
trial variable value pos group
1 1
Re: [R] comparing variances/distributions
Dear Ana This really depends on your scientific question. The two techniques you have shown do different things and there must be many more which could be applied. Michael On 17/06/2020 20:57, Ana Marija wrote: Hello, I have p values from two distributions, Pold and Pnew head(m) CHR POS MARKER Pnew Pold 1: 1 785989 rs2980300 0.1419 0.9521 2: 1 1130727 rs10907175 0.1022 0.4750 3: 1 1156131 rs2887286 0.3698 0.5289 4: 1 1158631 rs6603781 0.1929 0.2554 5: 1 1211292 rs6685064 0.6054 0.2954 6: 1 1478153 rs3766180 0.6511 0.5542 ... In order to compare those two distributions (QQ plots shown in attach) does it make sense to use: var.test(m$Pold, m$Pnew, alternative = "two.sided") F test to compare two variances data: m$Pold and m$Pnew F = 0.99937, num df = 1454159, denom df = 1454159, p-value = 0.7057 alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval: 0.9970808 1.0016750 sample estimates: ratio of variances 0.9993739 Or some other test makes more sense? Thanks Ana __ 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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] Convergence in Monte Carlo Simulation
Dear Edward
Every time you call your function powercrosssw() it resets the seed so
you must be calling it multiple times in some way.
Michael
On 14/06/2020 13:57, Phat Chau wrote:
Thank you Michael.
I will clarify some more. The function in the first part of the code that I
posted generates the simulated dataset for a cluster randomized trial from the
simstudy package.
I am not quite clear what you mean by placing it outside the loop. So the goal
here is to create n = 1000 independent datasets with different (randomly drawn
values from the specified normal distributions not shown) for all of the
parameters. What I have tried to do is place the seed at the very top of all my
code in the past, but what that does is it leads to the creation of a single
dataset that gets repeated over and over n = 1000 times. Hence, there ends up
being no variability in the data (and power estimates from the p-values given
the stated and required power).
Regarding the counter, is it correct in this instance that the loop will
continue until n = 1000 iterations have successfully converged? I am not so
concerned with counting failures.
Thank you.
Edward
On 2020-06-14, 6:46 AM, "Michael Dewey" wrote:
I am not 100% clear what your code is doing as it gets a bit wangled as
you posted in HTML but here are a couple of thoughts.
You need to set the seed outside any loops so it happens once and for all.
I would test after trycatch and keep a separate count of failures and
successes as the failure to converge must be meaningful about the
scientific question whatever that is. At the moment your count appears
to be in the correct place to count successes.
Michael
On 14/06/2020 02:50, Phat Chau wrote:
> Hello,
>
> I put together the following code and am curious about its correctness.
My first question relates to the Monte Carlo simulations – the goal is to continue
to iterate until I get n = 1000 simulations where the model successfully
converges. I am wondering if I coded it correctly below with the while loop. Is
the idea that the counter increments by one only if “model” does not return a
string?
>
> I would also like to know how I can create n = 1000 independent data
sets. I think to do this, I would have to set a random number seed via set.seed()
before the creation of each dataset. Where would I enter set.seed in the syntax
below? Would it be in the function (as indicated in red)?
>
> powercrosssw <- function(nclus, clsize) {
>
>set.seed()
>
>cohortsw <- genData(nclus, id = "cluster")
>cohortsw <- addColumns(clusterDef, cohortsw)
>cohortswTm <- addPeriods(cohortsw, nPeriods = 8, idvars = "cluster", perName =
"period")
>cohortstep <- trtStepWedge(cohortswTm, "cluster", nWaves = 4, lenWaves = 1,
startPer = 1, grpName = "Ijt")
>
>pat <- genCluster(cohortswTm, cLevelVar = "timeID", numIndsVar = clsize,
level1ID = "id")
>
>dx <- merge(pat[, .(cluster, period, id)], cohortstep, by = c("cluster",
"period"))
>dx <- addColumns(patError, dx)
>
>setkey(dx, id, cluster, period)
>
>dx <- addColumns(outDef, dx)
>
>return(dx)
>
> }
>
> i=1
>
> while (i < 1000) {
>
>dx <- powercrosssw()
>
>#Fit multi-level model to simulated dataset
>model5 <- tryCatch(lme(y ~ factor(period) + factor(Ijt), data = dx, random =
~1|cluster, method = "REML"),
> warning = function(w) { "warning" }
>)
>
>if (! is.character(model5)) {
>
> coeff <- coef(summary(model5))["factor(Ijt)1", "Value"]
> pvalue <- coef(summary(model5))["factor(Ijt)1", "p-value"]
> error <- coef(summary(model5))["factor(Ijt)1", "Std.Error"]
> bresult <- c(bresult, coeff)
> presult <- c(presult, pvalue)
> eresult <- c(eresult, error)
>
> i <- i + 1
>}
> }
>
> Thank you so much.
>
>
>
> [[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] Convergence in Monte Carlo Simulation
I am not 100% clear what your code is doing as it gets a bit wangled as
you posted in HTML but here are a couple of thoughts.
You need to set the seed outside any loops so it happens once and for all.
I would test after trycatch and keep a separate count of failures and
successes as the failure to converge must be meaningful about the
scientific question whatever that is. At the moment your count appears
to be in the correct place to count successes.
Michael
On 14/06/2020 02:50, Phat Chau wrote:
Hello,
I put together the following code and am curious about its correctness. My
first question relates to the Monte Carlo simulations – the goal is to continue
to iterate until I get n = 1000 simulations where the model successfully
converges. I am wondering if I coded it correctly below with the while loop. Is
the idea that the counter increments by one only if “model” does not return a
string?
I would also like to know how I can create n = 1000 independent data sets. I
think to do this, I would have to set a random number seed via set.seed()
before the creation of each dataset. Where would I enter set.seed in the syntax
below? Would it be in the function (as indicated in red)?
powercrosssw <- function(nclus, clsize) {
set.seed()
cohortsw <- genData(nclus, id = "cluster")
cohortsw <- addColumns(clusterDef, cohortsw)
cohortswTm <- addPeriods(cohortsw, nPeriods = 8, idvars = "cluster", perName =
"period")
cohortstep <- trtStepWedge(cohortswTm, "cluster", nWaves = 4, lenWaves = 1, startPer =
1, grpName = "Ijt")
pat <- genCluster(cohortswTm, cLevelVar = "timeID", numIndsVar = clsize, level1ID =
"id")
dx <- merge(pat[, .(cluster, period, id)], cohortstep, by = c("cluster",
"period"))
dx <- addColumns(patError, dx)
setkey(dx, id, cluster, period)
dx <- addColumns(outDef, dx)
return(dx)
}
i=1
while (i < 1000) {
dx <- powercrosssw()
#Fit multi-level model to simulated dataset
model5 <- tryCatch(lme(y ~ factor(period) + factor(Ijt), data = dx, random =
~1|cluster, method = "REML"),
warning = function(w) { "warning" }
)
if (! is.character(model5)) {
coeff <- coef(summary(model5))["factor(Ijt)1", "Value"]
pvalue <- coef(summary(model5))["factor(Ijt)1", "p-value"]
error <- coef(summary(model5))["factor(Ijt)1", "Std.Error"]
bresult <- c(bresult, coeff)
presult <- c(presult, pvalue)
eresult <- c(eresult, error)
i <- i + 1
}
}
Thank you so much.
[[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.
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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 in ordinal mixed effects
Dear Susana
Without your dat it is hard to say (and it would have helped to know
where mixor() comes from) but this almost always means that ne of your
parameters to the call is not what you thought it was so trying str(res)
might be enlightening. Also I do not see anywhere in your example where
you define na unless it is really NA and you did not copy it correctly.
Michael
On 09/06/2020 12:33, SUSANA ALBERICH MESA wrote:
Hi,
I'm trying to run an ordinal mixed effects model with Mixor command. I have 65
cases and repeated visits in 0, 6, 9, 12 and 18 months. My code is the
following:
cannabis<-c(datos$cannabis0, datos$cannabis6, datos$cannabis9,
datos$cannabis12, datos$cannabis18)
time<-c(rep(0, 65), rep(6, 65), rep(9, 65), rep(12, 65), rep(18, 65))
id<-c(rep(datos$id, 5))
group<-c(rep(datos$group, 5))
res<-data.frame(cbind(id, group, time, cannabis))
names(res)<-c("id", "group", "time", "cannabis")
res<-res[order(res$id),]
cannabismod<-mixor(cannabis~ time + as.factor(group), data=res, id=id, na.exclude,
which.random.slope=na, link="logit")
summary(cannabismod)
However, I have obtained this error:
Error in xj[i] : invalid subscript type 'closure'
Please, could anyone help me to solve it?
Many thanks,
Susana
[https://edukiak.osakidetza.net/coronavirus/pie_email7.jpg]
[[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.
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] a question of etiquette
You might get better answers on the list dedicated to package development r-pkg-devel This may have already been discussed there so a quick look at the archive might also help you. On 01/06/2020 17:34, Adelchi Azzalini wrote: The new version of a package which I maintain will include a new function which I have ported to R from Matlab. The documentation of this R function indicates the authors of the original Matlab code, reference to their paper, URL of the source code. Question: is this adequate, or should I include them as co-authors of the package, or as contributors, or what else? Is there a general policy about this matter? Adelchi Azzalini http://azzalini.stat.unipd.it/ __ 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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] struggling with apply
This is like "Name that Tune." Can anyone do it in FEWER characters? :-) On May 27, 2020, at 4:32 PM, Mathew Guilfoyle wrote: A bit quicker: t(pmin(t(somematrix), UB)) On 27 May 2020, at 20:56, Bert Gunter mailto:bgunter.4...@gmail.com>> wrote: Jeff: Check it! somematrix <- matrix(c(1,4,3,6,3,9,12,8,5,7,11,11),nrow=3,ncol=4) UB=c(2.5, 5.5, 8.5, 10.5) apply( somematrix, 2, function( x ) pmin( x, UB ) ) [,1] [,2] [,3] [,4] [1,]1 2.5 2.5 2.5 [2,]4 3.0 5.5 5.5 [3,]3 8.5 5.0 8.5 [4,]1 6.0 10.5 7.0 Not what was wanted. Am I missing something? 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 Wed, May 27, 2020 at 12:38 PM Jeff Newmiller mailto:jdnew...@dcn.davis.ca.us>> wrote: Sigh. Transpose? apply( somematrix, 2, function( x ) pmin( x, UB ) ) On May 27, 2020 11:22:06 AM PDT, Bert Gunter mailto:bgunter.4...@gmail.com>> wrote: Better, I think (no indexing): t(apply(somematrix,1,function(x)pmin(x,UB))) 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 Wed, May 27, 2020 at 10:56 AM Rui Barradas mailto:ruipbarra...@sapo.pt>> wrote: Hello, Try pmin. And loop by column/UB index with sapply/seq_along. sapply(seq_along(UB), function(i) pmin(UB[i], somematrix[,i])) # [,1] [,2] [,3] [,4] #[1,] 1.0 5.5 8.5 7.0 #[2,] 2.5 3.0 8.0 10.5 #[3,] 2.5 5.5 5.0 10.5 Hope this helps, Rui Barradas Às 18:46 de 27/05/20, Michael Ashton escreveu: Hi - I have a matrix of n rows and 4 columns. I want to cap the value in each column by a different upper bound. So, suppose my matrix is somematrix <- matrix(c(1,4,3,6,3,9,12,8,5,7,11,11),nrow=3,ncol=4) somematrix [,1] [,2] [,3] [,4] [1,]16 127 [2,]438 11 [3,]395 11 Now I want to have the maximum value in each column described by UB=c(2.5, 5.5, 8.5, 10.5) So that the right answer will look like: [,1] [,2][,3] [,4] [1,]1 5.5 8.57 [2,]2.538 10.5 [3,]2.5 5.5 510.5 I've tried a few things, like: newmatrix <- apply(somematrix,c(1,2),function(x) min(UB,x)) but I can't figure out to apply the relevant element of the UB list to the right element of the matrix. When I run the above, for example, it takes min(UB,x) over all UB, so I get: newmatrix [,1] [,2] [,3] [,4] [1,] 1.0 2.5 2.5 2.5 [2,] 2.5 2.5 2.5 2.5 [3,] 2.5 2.5 2.5 2.5 I'm sure there's a simple and elegant solution but I don't know what it is! Thanks in advance, Mike Michael Ashton, CFA Managing Principal Enduring Investments LLC W: 973.457.4602 C: 551.655.8006 Schedule a Call: https://calendly.com/m-ashton [[alternative HTML version deleted]] __ R-help@r-project.org<mailto: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. [[alternative HTML version deleted]] __ R-help@r-project.org<mailto: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. [[alternative HTML version deleted]] __ R-help@r-project.org<mailto: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<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 -- 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] struggling with apply
Always amazes me how many ways there are to do these things, none of which I was able to find myself. Thanks! I think the key here was ‘pmin,’ which I didn’t know before. Michael Ashton, CFA Managing Principal Enduring Investments LLC W: 973.457.4602 C: 551.655.8006 Schedule a Call: https://calendly.com/m-ashton From: Bert Gunter [mailto:bgunter.4...@gmail.com] Sent: Wednesday, May 27, 2020 2:22 PM To: Rui Barradas Cc: Michael Ashton; r-help@r-project.org Subject: Re: [R] struggling with apply Better, I think (no indexing): t(apply(somematrix,1,function(x)pmin(x,UB))) 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 Wed, May 27, 2020 at 10:56 AM Rui Barradas mailto:ruipbarra...@sapo.pt>> wrote: Hello, Try pmin. And loop by column/UB index with sapply/seq_along. sapply(seq_along(UB), function(i) pmin(UB[i], somematrix[,i])) # [,1] [,2] [,3] [,4] #[1,] 1.0 5.5 8.5 7.0 #[2,] 2.5 3.0 8.0 10.5 #[3,] 2.5 5.5 5.0 10.5 Hope this helps, Rui Barradas Às 18:46 de 27/05/20, Michael Ashton escreveu: > Hi - > > I have a matrix of n rows and 4 columns. > > I want to cap the value in each column by a different upper bound. So, > suppose my matrix is > > somematrix <- matrix(c(1,4,3,6,3,9,12,8,5,7,11,11),nrow=3,ncol=4) >> somematrix > [,1] [,2] [,3] [,4] > [1,]16 127 > [2,]438 11 > [3,]395 11 > > Now I want to have the maximum value in each column described by > UB=c(2.5, 5.5, 8.5, 10.5) > > So that the right answer will look like: > [,1] [,2][,3] [,4] > [1,]1 5.5 8.57 > [2,]2.538 10.5 > [3,]2.5 5.5 510.5 > > I've tried a few things, like: > newmatrix <- apply(somematrix,c(1,2),function(x) min(UB,x)) > > but I can't figure out to apply the relevant element of the UB list to the > right element of the matrix. When I run the above, for example, it takes > min(UB,x) over all UB, so I get: > > newmatrix > [,1] [,2] [,3] [,4] > [1,] 1.0 2.5 2.5 2.5 > [2,] 2.5 2.5 2.5 2.5 > [3,] 2.5 2.5 2.5 2.5 > > I'm sure there's a simple and elegant solution but I don't know what it is! > > Thanks in advance, > > Mike > > Michael Ashton, CFA > Managing Principal > > Enduring Investments LLC > W: 973.457.4602 > C: 551.655.8006 > Schedule a Call: https://calendly.com/m-ashton > > > [[alternative HTML version deleted]] > > __ > R-help@r-project.org<mailto: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<mailto: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. [[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] struggling with apply
That's it! Thanks. Learn something new every day! Michael Ashton, CFA Managing Principal Enduring Investments LLC W: 973.457.4602 C: 551.655.8006 Schedule a Call: https://calendly.com/m-ashton -Original Message- From: Rui Barradas [mailto:ruipbarra...@sapo.pt] Sent: Wednesday, May 27, 2020 1:51 PM To: Michael Ashton; r-help@r-project.org Subject: Re: [R] struggling with apply Hello, Try pmin. And loop by column/UB index with sapply/seq_along. sapply(seq_along(UB), function(i) pmin(UB[i], somematrix[,i])) # [,1] [,2] [,3] [,4] #[1,] 1.0 5.5 8.5 7.0 #[2,] 2.5 3.0 8.0 10.5 #[3,] 2.5 5.5 5.0 10.5 Hope this helps, Rui Barradas Às 18:46 de 27/05/20, Michael Ashton escreveu: > Hi - > > I have a matrix of n rows and 4 columns. > > I want to cap the value in each column by a different upper bound. So, > suppose my matrix is > > somematrix <- matrix(c(1,4,3,6,3,9,12,8,5,7,11,11),nrow=3,ncol=4) >> somematrix > [,1] [,2] [,3] [,4] > [1,]16 127 > [2,]438 11 > [3,]395 11 > > Now I want to have the maximum value in each column described by > UB=c(2.5, 5.5, 8.5, 10.5) > > So that the right answer will look like: > [,1] [,2][,3] [,4] > [1,]1 5.5 8.57 > [2,]2.538 10.5 > [3,]2.5 5.5 510.5 > > I've tried a few things, like: > newmatrix <- apply(somematrix,c(1,2),function(x) min(UB,x)) > > but I can't figure out to apply the relevant element of the UB list to the > right element of the matrix. When I run the above, for example, it takes > min(UB,x) over all UB, so I get: > > newmatrix > [,1] [,2] [,3] [,4] > [1,] 1.0 2.5 2.5 2.5 > [2,] 2.5 2.5 2.5 2.5 > [3,] 2.5 2.5 2.5 2.5 > > I'm sure there's a simple and elegant solution but I don't know what it is! > > Thanks in advance, > > Mike > > Michael Ashton, CFA > Managing Principal > > Enduring Investments LLC > W: 973.457.4602 > C: 551.655.8006 > Schedule a Call: https://calendly.com/m-ashton > > > [[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.
[R] struggling with apply
Hi - I have a matrix of n rows and 4 columns. I want to cap the value in each column by a different upper bound. So, suppose my matrix is somematrix <- matrix(c(1,4,3,6,3,9,12,8,5,7,11,11),nrow=3,ncol=4) > somematrix [,1] [,2] [,3] [,4] [1,]16 127 [2,]438 11 [3,]395 11 Now I want to have the maximum value in each column described by UB=c(2.5, 5.5, 8.5, 10.5) So that the right answer will look like: [,1] [,2][,3] [,4] [1,]1 5.5 8.57 [2,]2.538 10.5 [3,]2.5 5.5 510.5 I've tried a few things, like: newmatrix <- apply(somematrix,c(1,2),function(x) min(UB,x)) but I can't figure out to apply the relevant element of the UB list to the right element of the matrix. When I run the above, for example, it takes min(UB,x) over all UB, so I get: newmatrix [,1] [,2] [,3] [,4] [1,] 1.0 2.5 2.5 2.5 [2,] 2.5 2.5 2.5 2.5 [3,] 2.5 2.5 2.5 2.5 I'm sure there's a simple and elegant solution but I don't know what it is! Thanks in advance, Mike Michael Ashton, CFA Managing Principal Enduring Investments LLC W: 973.457.4602 C: 551.655.8006 Schedule a Call: https://calendly.com/m-ashton [[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] text annotation on Manhattn plot in qqman
Dear Ana That looks like something is hard coded in manhattan(). The simplest thing might be to contact the maintainer of the package and ask. You can make a copy of manhattan or textxy and modify them but I think the maintainer is the simplest course of action. Michael On 20/05/2020 16:10, Ana Marija wrote: HI Michael, Thank you so much! That worked!!! Now I am just trying to increase the size of text of SNP and GENE on plot I tried this: a$newname <- paste(a$SNP,"\n", a$GENE) manhattan(a, chr="CHR", bp="BP", snp="newname", p="P",cex = 0.5,annotatePval = 0.0001) but I am getting this error: Error in textxy(topSNPs$pos, -log10(topSNPs$P), offset = 0.625, labs = topSNPs$SNP, : formal argument "cex" matched by multiple actual arguments Do you by any chance know how to do this? Cheers Ana On Wed, May 20, 2020 at 4:12 AM Michael Dewey wrote: a$newname <- paste(a$SNP, a$GENE) manhattan(a, chr="CHR", bp="BP", snp="newname", p="P",annotatePval = 0.0001) However note that I do not use manhattan() so you may need to alter the parameters as I am assuming a is where it finds the remaining parameters. You may also need to play with the sep =, and collapse = parameters to paste() to get the precise layout you want. Michael On 19/05/2020 17:21, Ana Marija wrote: Hi Michael, can you please send me code how that would be done? Thanks Ana On Tue, May 19, 2020 at 11:18 AM Michael Dewey wrote: Dear Ana Perhaps paste together SNP and GENE using paste() and then supply that as the snp parameter. Michael On 19/05/2020 17:12, Ana Marija wrote: Hello, I am making manhattan plot with: library(qqman) manhattan(a, chr="CHR", bp="BP", snp="SNP", p="P",annotatePval = 0.0001) and I would like to annotate these two SNPs which are above the threshold so that they have GENE name beside them: a[a$SNP=="rs4081570",] SNPP CHR BP GENE 1 rs4081570 6.564447e-05 19 15918647 UCA1 a[a$SNP=="rs11867934",] SNPP CHR BP GENE 1021 rs11867934 6.738066e-06 17 16933404 FLCN Right now my plot only has SNP name for those 2, how to add GENE names (FLCN and UCA1 as well) Please advise Ana __ 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. -- Michael http://www.dewey.myzen.co.uk/home.html -- Michael http://www.dewey.myzen.co.uk/home.html -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] text annotation on Manhattn plot in qqman
a$newname <- paste(a$SNP, a$GENE) manhattan(a, chr="CHR", bp="BP", snp="newname", p="P",annotatePval = 0.0001) However note that I do not use manhattan() so you may need to alter the parameters as I am assuming a is where it finds the remaining parameters. You may also need to play with the sep =, and collapse = parameters to paste() to get the precise layout you want. Michael On 19/05/2020 17:21, Ana Marija wrote: Hi Michael, can you please send me code how that would be done? Thanks Ana On Tue, May 19, 2020 at 11:18 AM Michael Dewey wrote: Dear Ana Perhaps paste together SNP and GENE using paste() and then supply that as the snp parameter. Michael On 19/05/2020 17:12, Ana Marija wrote: Hello, I am making manhattan plot with: library(qqman) manhattan(a, chr="CHR", bp="BP", snp="SNP", p="P",annotatePval = 0.0001) and I would like to annotate these two SNPs which are above the threshold so that they have GENE name beside them: a[a$SNP=="rs4081570",] SNPP CHR BP GENE 1 rs4081570 6.564447e-05 19 15918647 UCA1 a[a$SNP=="rs11867934",] SNPP CHR BP GENE 1021 rs11867934 6.738066e-06 17 16933404 FLCN Right now my plot only has SNP name for those 2, how to add GENE names (FLCN and UCA1 as well) Please advise Ana __ 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. -- Michael http://www.dewey.myzen.co.uk/home.html -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] text annotation on Manhattn plot in qqman
Dear Ana Perhaps paste together SNP and GENE using paste() and then supply that as the snp parameter. Michael On 19/05/2020 17:12, Ana Marija wrote: Hello, I am making manhattan plot with: library(qqman) manhattan(a, chr="CHR", bp="BP", snp="SNP", p="P",annotatePval = 0.0001) and I would like to annotate these two SNPs which are above the threshold so that they have GENE name beside them: a[a$SNP=="rs4081570",] SNPP CHR BP GENE 1 rs4081570 6.564447e-05 19 15918647 UCA1 a[a$SNP=="rs11867934",] SNPP CHR BP GENE 1021 rs11867934 6.738066e-06 17 16933404 FLCN Right now my plot only has SNP name for those 2, how to add GENE names (FLCN and UCA1 as well) Please advise Ana __ 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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] the volcano orientation
Does anyone know why 'volcano' is oriented as it is?
image(volcano) ## filled.contour is the same
I know it's all arbitrary, but north-up is a pretty solid convention. Is
there any reason why the classic 'image()' example data set would not
default to this orientation?
A Google map of the site (in Web Mercator):
https://www.google.com/maps/place/Maungawhau+%2F+Mount+Eden/@-36.8763271,174.7619561,856m/data=!3m1!1e3!4m8!1m2!2m1!1smaungawhau!3m4!1s0x6d0d47db8d7bd1ff:0x8bcffe2a5c7360d2!8m2!3d-36.867!4d174.767
For image(), the north-up orientation is 't(volcano[,ncol(volcano):1])'.
If you are interested in a roughly georeferenced version I have code here:
https://gist.github.com/mdsumner/20fe3ffa04421bf8e0517c19085e5fd8
(Also see fortunes::fortune("conventions") )
Best, Mike
--
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsum...@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] NA command in a 'for' loop
Just a thought Helen but is x being treated as a real and what you think
are zero and are printed as zero are in fact some very small number? If
so you need to alter your test appropriately.
Michael
On 20/04/2020 17:25, Helen Sawaya wrote:
Thank you for your reply.
I tried d[] <- lapply(d, function(x) {is.na(x) <- x == 0; x})
but I am still getting zeros instead of NAs in my output..
I wonder if the problem is that some of my data files don't have any zeros
(participants made no errors)..
From: Rui Barradas
Sent: Monday, April 20, 2020 9:05 AM
To: Helen Sawaya ; r-help@R-project.org
Subject: Re: [R] NA command in a 'for' loop
Hello,
Instead of
d[d == 0] <- NA
try
d[] <- lapply(d, function(x) {is.na(x) <- x == 0; x})
Also, in the first for loop
paste(i, sep = "")
does nothing, it's the same as i.
And the same for
(d2$V4 == 1) == TRUE
Since (d2$V4 == 1) already is FALSE/TRUE there is no need for
(.) == TRUE
Hope this helps,
Rui Barradas
Às 20:52 de 19/04/20, Helen Sawaya escreveu:
Dear R experts,
I am using a 'for' loop to apply commands to multiple datasets (each file is one
participant). The only one not working is the command that identifies zeros in my
datasets and changes them to NAs. But when I look at the output, zeros ("0")
are still present. Surprisingly, the functions work fine when I apply them to a single
dataset (outside the loop). I've tried:
all.files <- list.files(".")
txt.files <- grep("threat.txt",all.files,value=T)
for(i in txt.files){
d <- read.table(paste(i,sep=""),header=F)
d[d==0] <- NA #replace zeros with NA
write.table(d, paste0(i,".tlbs.txt"), quote=FALSE, row.names=TRUE)}
d<-d[ ,-c(10,11)]
d2<-d[complete.cases(d), ]
d2$V4<-as.numeric(d2$V4)
congruent <- (d2$V4 == 1) == TRUE
x <- get_tlbs(d2$V14, congruent, prior_weights = NULL, method = "weighted",
fill_gaps = FALSE)
write.table(x, paste0(i,".tlbs.txt"), quote=FALSE, row.names=TRUE)}
I've also tried:
for(i in txt.files){
d <- read.table(paste(i,sep=""),header=F)
if (0 %in% d)
{replace_with_na(d,replace = list(x = 0))} # replace zeros with NA
d<-d[ ,-c(10,11)]
d2<-d[complete.cases(d), ]
d2$V4<-as.numeric(d2$V4)
congruent <- (d2$V4 == 1) == TRUE
x <- get_tlbs(d2$V14, congruent, prior_weights = NULL, method = "weighted",
fill_gaps = FALSE)
write.table(x, paste0(i,".summaryoutput.txt"), quote=FALSE, row.names=TRUE)}
Thank you for your help.
Sincerely
Helen
[[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.
[[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.
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] how to get all strings in a data frame that start with a particular string
Dear Ana
Would it not be possible to use grep instead of grepl and get the values
using the value = TRUE parameter?
Michael
On 10/04/2020 17:15, Ana Marija wrote:
Hello,
Hello,
I have a data frame (tot) with about 2000 columns. How can I extract
from it all strings that start with E14?
I tried this:
e14 <- sapply(tot, function(x) grepl("^E14", x))
but this returns me just TRUE and FALSE vector, how do I get actual
strings that start with E14?
Thanks
Ana
__
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.
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] Same results but different functions ?
The documentation suggests that the rlm method for a formula does not
have psi as a parameter. Perhaps try using the method for a matrix x and
a vector y.
Michael
On 23/03/2020 12:39, varin sacha via R-help wrote:
Dear R-experts,
The rlm command in the MASS package command implements several versions of
robust regression, for example the Huber and the Tukey (bisquare weighting
function) estimators.
In my R code here below I try to get the Tukey (bisquare weighting function)
estimation, R gives me an error message : Error in statistic(data, original,
...) : unused argument (psi = psi.bisquare)
If I cancel psi=psi.bisquare my code is working but IMHO I will get the Huber
estimation and not the Tukey. So how can I get the Tukey ? Many thanks for your
help.
# # # # # # # # # # # # # # # # # # # # # # # #
install.packages( "boot",dependencies=TRUE )
install.packages( "MASS",dependencies=TRUE )
library(boot)
library(MASS)
n<-50
b<-runif(n, 0, 5)
z <- rnorm(n, 2, 3)
a <- runif(n, 0, 5)
y_model<- 0.1*b - 0.5 * z - a + 10
y_obs <- y_model +c( rnorm(n*0.9, 0, 0.1), rnorm(n*0.1, 0, 0.5) )
df<-data.frame(b,z,a,y_obs)
# function to obtain MSE
MSE <- function(data, indices, formula) {
d <- data[indices, ] # allows boot to select sample
fit <- rlm(formula, data = d)
ypred <- predict(fit)
d[["y_obs "]] <-y_obs
mean((d[["y_obs"]]-ypred)^2)
}
# Make the results reproducible
set.seed(1234)
# bootstrapping with 600 replications
results <- boot(data = df, statistic = MSE,
R = 600, formula = y_obs ~ b+z+a, psi = psi.bisquare)
str(results)
boot.ci(results, type="bca" )
# # # # # # # # # # # # # # # # # # # # # # # # #
__
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.
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] how do I add text lables on QQ plot
Dear Ana
You can specify the first three parameters to text() as vectors so it is
all done in one call. That may or may not answer your question.
Michael
On 12/03/2020 14:08, Ana Marija wrote:
HI David,
thank you for getting back to me.
Is there is a way for qq() to pick up text label names on its own or I
have to specify each one manually?
like in this example:
text( 2, 6, "arbitrary")
this is dput for
a=head(fdr2_sorted)
dput(a)
structure(list(NAME = c("GO_DNA_PACKAGING_COMPLEX", "GO_PROTEIN_DNA_COMPLEX",
"GO_RESPONSE_TO_TYPE_I_INTERFERON", "GO_RESPONSE_TO_INTERFERON_GAMMA",
"GO_CELLULAR_RESPONSE_TO_INTERFERON_GAMMA", "GO_GRANULOCYTE_MIGRATION"
), `GS follow link to MSigDB` = c("GO_DNA_PACKAGING_COMPLEX",
"GO_PROTEIN_DNA_COMPLEX", "GO_RESPONSE_TO_TYPE_I_INTERFERON",
"GO_RESPONSE_TO_INTERFERON_GAMMA", "GO_CELLULAR_RESPONSE_TO_INTERFERON_GAMMA",
"GO_GRANULOCYTE_MIGRATION"), `GS DETAILS` = c("Details ...",
"Details ...", "Details ...", "Details ...", "Details ...", "Details ..."
), SIZE = c(77L, 132L, 52L, 101L, 85L, 43L), ES = c(0.6757226,
0.59581786, -0.7521569, -0.63704145, -0.6571892, -0.7332099),
NES = c(2.466745, 2.3465202, -2.5331483, -2.4204733, -2.4021535,
-2.3989832), `NOM p-val` = c(0, 0, 0, 0, 0, 0), `FDR q-val` = c(0,
0, 0, 0, 0, 0), `FWER p-val` = c(0, 0, 0, 0, 0, 0), `RANK AT MAX` =
c(L,
1516L, 1427L, 1819L, 1216L, 491L), `LEADING EDGE` = c("tags=43%,
list=10%, signal=47%",
"tags=39%, list=13%, signal=45%", "tags=54%, list=12%, signal=61%",
"tags=45%, list=16%, signal=52%", "tags=38%, list=11%, signal=42%",
"tags=28%, list=4%, signal=29%"), V12 = c(NA, NA, NA, NA,
NA, NA), FDR.q.val2 = c(1e-10, 1e-10, 1e-10, 1e-10, 1e-10,
1e-10)), class = c("data.table", "data.frame"), row.names = c(NA,
-6L), .internal.selfref = )
b=head(fdr1_sorted)
dput(b)
structure(list(NAME = c("GO_DNA_PACKAGING_COMPLEX", "GO_PROTEIN_DNA_COMPLEX",
"GO_RESPONSE_TO_TYPE_I_INTERFERON", "GO_RESPONSE_TO_INTERFERON_GAMMA",
"GO_CELLULAR_RESPONSE_TO_INTERFERON_GAMMA", "GO_GRANULOCYTE_MIGRATION"
), `GS follow link to MSigDB` = c("GO_DNA_PACKAGING_COMPLEX",
"GO_PROTEIN_DNA_COMPLEX", "GO_RESPONSE_TO_TYPE_I_INTERFERON",
"GO_RESPONSE_TO_INTERFERON_GAMMA", "GO_CELLULAR_RESPONSE_TO_INTERFERON_GAMMA",
"GO_GRANULOCYTE_MIGRATION"), `GS DETAILS` = c("Details ...",
"Details ...", "Details ...", "Details ...", "Details ...", "Details ..."
), SIZE = c(77L, 132L, 52L, 101L, 85L, 43L), ES = c(0.6757226,
0.59581786, -0.7521569, -0.63704145, -0.6571892, -0.7332099),
NES = c(2.466745, 2.3465202, -2.5331483, -2.4204733, -2.4021535,
-2.3989832), `NOM p-val` = c(0, 0, 0, 0, 0, 0), `FDR q-val` = c(0,
0, 0, 0, 0, 0), `FWER p-val` = c(0, 0, 0, 0, 0, 0), `RANK AT MAX` =
c(L,
1516L, 1427L, 1819L, 1216L, 491L), `LEADING EDGE` = c("tags=43%,
list=10%, signal=47%",
"tags=39%, list=13%, signal=45%", "tags=54%, list=12%, signal=61%",
"tags=45%, list=16%, signal=52%", "tags=38%, list=11%, signal=42%",
"tags=28%, list=4%, signal=29%"), V12 = c(NA, NA, NA, NA,
NA, NA), group = c(2, 2, 4, 4, 4, 4), FDR.q.val2 = c(1e-10,
1e-10, 1e-10, 1e-10, 1e-10, 1e-10)), class = c("data.table",
"data.frame"), row.names = c(NA, -6L), .internal.selfref = )
library(qqman)
qq(fdr2_sorted$FDR.q.val2, main = "RG_All", pch = 16,
col=fdr1_sorted$group, cex = 0.8, las = 1)
Please advise
On Wed, Mar 11, 2020 at 11:21 PM David Winsemius wrote:
On 3/10/20 9:51 PM, Ana Marija wrote:
Hello,
I am making QQ plot via:
library(ggman)
qq(fdr2_sorted$FDR.q.val2, main = "RG_All", pch = 17,
col=fdr1_sorted$group, cex = 1, las = 1)
I think you may be confusing the audience. There is no qq function in
the ggman package. There is however a qq function in the qqman package.
Running the example in help page for qqman::qq and looking at the code
suggests this is a base plot function, so the text function will allow
you to put any particular string within the plot area:
library(qqman)
qq(gwasResults$P)
text( 2, 6, "arbitrary") # puts text "arbitrary" at postion (x=2, y=6)
data frames used look like this:
head(fdr1_sorted)
You should use `dput` to post reproducible data examples.
HTH;
David.
NAME GS follow
link to MSigDB GS DETAILS SIZE ES NES NOM p-val
1: GO_DNA_PACKAGING_COMPLEX
GO_DNA_PACKAGING_COMPLEX Details ...
[R] Announcement: Summer Statistics Institute at UT Austin
The Department of Statistics and Data Sciences at The University of Texas at Austin is hosting the 13th annual UT Summer Statistics Institute (SSI) May 26–29, 2020. SSI OVERVIEW: * SSI attracts participants from academia, health professions and marketing firms. * Participants acquire new statistical knowledge and skills through hands-on data analysis. * SSI provides an exciting venue to build statistical knowledge alongside a diverse and dynamic audience. UT's Summer Statistics Institute (SSI) offers intensive four-day workshops on diverse topics from introductory data sciences to advanced statistics. Whether you are new to data analysis or a seasoned statistician, SSI provides a unique hands-on opportunity to acquire valuable skills directly from experts in the field. The UT Summer Statistics Institute (SSI) is open to 700 participants. Courses span introductory statistics, statistical software, statistical methods and statistics applications. Each course will meet for four half-days, either mornings or afternoons, for a total of twelve hours. Instructors will post lectures, datasets, exercises, and course information on a website accessible to enrolled participants. There will be no examinations, and participants will receive certificates upon completion. Academic credit will not be issued. Please carefully check the specified prerequisite knowledge before enrolling in a course. R Specific Courses for this year include: * Introduction to Data Analysis and Graphics Using R (AM) * Introduction to Data Analysis and Graphics Using R (PM) * Informed Decision Making from Data: Regression Analysis * Introduction to Bayesian Statistics The Department of Statistics and Data Sciences now offers credit for educational professional development opportunities during our 2020 Summer Statistics Institute! Teachers currently working in PK–12 settings can earn 12 hours of Continuing Professional Education (CPE) by attending one of our 25 exciting SSI courses held May 26-29, 2020. Courses will be held on the UT Campus in Patton Hall (RLP), Flawn Academic Center (FAC) and Robert A. Welch Hall (WEL). For more information and to register: https://stat.utexas.edu/training/ssi - MICHAEL J. MAHOMETA Director of Statistical Consulting and Professional Education Department of Statistics and Data Sciences | The University of Texas at Austin 2317 Speedway, Stop D9800 | Austin, TX 78712 o: 512-471-4542 | stat.utexas.edu __ 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] R CMD check report error not finding function
Dear Ruben
I do not think your export pattern matches .foo so that may be the problem.
Michael
On 05/03/2020 09:10, Ruben wrote:
Dear R-help listers,
I am creating a package for a friend using his scripts and I run into a
problem when doing the check.
Inside a macro created with gtools::defmacro(par1, expr= ... there is a
call to stats::optim and the function for the optimization is called
.foo(), which is listed as one of the functions to be included in the
package, not to be documented, using package.skeleton. R CMD check
reports that it cannot find .foo(). This runs OK as a script but not as
a package. I believe it has to do with namespaces. The NAMESPACE file is
just this:
exportPattern("^[[:alpha:]]+")
importFrom("stats", "optim")
import(gtools)
importFrom("gtools", "defmacro")
Any help would be greatly appreciated
Ruben
--
Ruben H. Roa-Ureta, Ph. D.
Consultant, ORCID ID -0002-9620-5224
Center for Environment and Water, Research Institute,
King Fahd University of Petroleum and Minerals,
KFUPM Box 1927, Dhahran 31261, Saudi Arabia
Office Phone : 966-3-860-7850
Cellular Phone : 966-540026401
إن المعلومات الواردة في هذا البريد الإلكتروني ومرفقاته (إن وجدت)، قد
تكون خاصة أو سرية؛ فإذا لم تكن المقصود بهذه الرسالة؛ فيُرجى منك حذفها
ومرفقاتها من نظامك وإخطار المرسل بخطأ وصولها إليك فورا. كما لا يجوز نسخ
أي جزء منها أو مرفقاتها ، أو الإفصاح عن محتوياتها لأي شخص أو استعمالها
لأي غرض آخر. إن جامعة الملك فهد للبترول والمعادن لا تتحمل مسؤولية
التغييرات التي يتم إجراؤها على هذه الرسالة بعد إرسالها. وإن البيانات أو
الآراء المعبر عنها في هذا البريد، هي بيانات تخص مُرسلها، ولا تعكس
بالضرورة رأي وبيانات الجامعة. كما لا تتحمل الجامعة مسؤولية أي تأثير ينتج
عن هذه الرسالة أوعن أي فيروس قد تحمله.
DISCLAIMER: The information in this email and its attach...{{dropped:11}}
__
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.
--
Michael
http://www.dewey.myzen.co.uk/home.html
__
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] testing my package : unstated dependency to self in package tests
When something similar happened to me I found it went away when I added Suggests: to the DESCRIPTION file. Whether this will work for you I have no idea. Michael On 16/02/2020 11:03, Servet Ahmet Çizmeli wrote: I am updating my CRAN package geoSpectral. I get the following Warning during R CMD check : ... * checking for unstated dependencies in �tests� ... WARNING 'library' or 'require' call not declared from: �geoSpectral� All the .R files I have under the testhat directory begin by : library(geoSpectral) library(testthat) and there I call package functions directly (without the prefix geoSpectal:: ) See https://github.com/cran/geoSpectral/blob/master/tests/testthat/Spectra_tests.R Searching the web, I found examples where the same Warning has been issued for some other packages. But in my case the package in question is my own package I am testing Confused and at loss. Anyone with ideas? regards Servet [[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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] Information about pairwise Wilcoxon Sum Rank Tests
Dear Gabriela Apologies if you have already tried this but: 1 - you can see the code which it uses by typing pairwise.wilcox.test at the command line. 2 - to find you more about the method of adjustment you need to follow the documentation for p.adjust 3 - you may also need to look into wilcox.test which pairwise.wilcox.test uses. If that does not help come back and tell us more about how far you have got so people know what else you need. Michael On 06/02/2020 10:04, Gabriela Hoff wrote: Dear Sir or Madam My name is Gabriela Hoff and I and an R user. I am facing problems in finding detailed information about pairwise Wilcoxon Rank Sum Tests (pairwise.wilcox.test). I would like to understand better the way this test is performed and the calculation of p-value, especially the procedure used on the corrections for multiple testing. I do appreciate it a lot if you could send me information about that since I was unable to find those details in the R documentation. Best regards, -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] R package for meta-analysis from z scores
Dear James Your question really boil down to whether you can estimate tau^2, the between study variance of the effect sizes, if you only have p-values. As far as I can see the answer has to be no. Michael On 16/01/2020 13:10, james poweraid wrote: Hello, I have a set of Z scores from an N=2 studies and I need to run a meta-analysis across the two Z scores for many N variables. I do not have effect sizes and SEs. I realize there are many different meta analysis packages in R, but I only have Z scores and it seems to me this is limiting. I am using the metap package because this very conveniently accepts directly p values which I can get from my Z scores. I have applied what is called in metap package the sumlog Fisher’s method Chi square (2 df). I have three questions: 1. Is this the same as a fixed effect meta-analysis without weights? 2. Is there any way to do a random effects meta-analysis starting only with Z scores? 3. Is there a way to get the I2 heterogeneity across studies from this or from other packages? It looks like I need the estimates for this. Is there another package that can easily get this using as input Z scores or P values? Data library(metap) zscores = cbind(c(4, 2, 0.1),c(5, 2.5, 0.1)) pvalues = apply(zscores, 1:2, function(x) 2*(pnorm( abs(x) , lower.tail=F))) meta = apply(pvalues, MARGIN=1, FUN=sumlog) Thank you much [[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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] "In sqrt(VS) : NaNs produced"
Your script and data were stripped so we are none the wiser I am afraid. You need to embed the script in your post and give a minimal data-set which exhibits the problem using dput() and embed that in the post too. Michael On 19/01/2020 08:25, Atul Saini wrote: Hello R, I am attaching the script and data please help me to solve the problem of "In sqrt(VS) : NaNs produced" with the p value of dumy$Mar Regards, __ 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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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] all the MAE metric values are missing (Error message)
What Jim is alluding to is that sometimes in the process of reading in data a small typo can mean that what was intended to be a numeric variable is read in as a factor. So he was suggesting that you double check that this has not happened to you. Michael On 23/12/2019 11:45, Neha gupta wrote: Hi Jim, Another possibility is that the values used in the initial calculation have been read in as factors Which calculation you are talking about? I did not use factors as variable. Regards On Mon, Dec 23, 2019 at 3:12 AM Jim Lemon wrote: Hi Neha, Well, that's a clue to why you are getting NAs: log10(0) [1] -Inf Another possibility is that the values used in the initial calculation have been read in as factors. Jim On Mon, Dec 23, 2019 at 10:55 AM Neha gupta wrote: Hi Jim The objective function is passed to san_res where we have defined the 4 parameters of gbm and the values are initialized in san_res. The output variable price has only three values: 0, 1, 2 (like categorical values), so someone told me try to remove the log10 from the price. [[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. -- Michael http://www.dewey.myzen.co.uk/home.html __ 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.
