Re: [R] Displaying Compositional Data With 46 Parts
On Fri, 17 Jul 2015, Aaron Mackey wrote: One immediate question is how independent you believe the 46 components to be, and whether certain components could be reduced or otherwise coordinately-modeled; a heatmap of your 46x46 pairwise correlations should be informative. Also consider log-scaling, especially if some component fractions can be very small/minor components compared to large/major components. Aaron, Most components are elements spread across the periodic table; the rest are composits such as total dissolved/suspended solids, acid neutralizing capacity, specific conductance, bicarbonate. A heat map will be drawn. Log scaling is the key to working with compositional data. The log-ratios can be calculated using three equations with two being most frequently applied, each depending on the analysis to follow. Thanks, Rich __ 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] Displaying Compositional Data With 46 Parts
On Fri, 17 Jul 2015, John Kane wrote: I think this is more a technical question for the subject matter experts than for R-help if I am understanding the question correctly. John, I agree completely. Unfortunately, there is no R SIG devoted to CoDA, nor any other mail list or Web forum that I've been able to find. There is a CoDA Web site but no active forum. That said, there seems to be an R package called compositional and a corresponding book http://www.springer.com/us/book/9783642368080 that may help. If nothing else the names of the author(s) of the package or book or the references may give you some leads on some decent papers. In addition to compositions, there's robCompositions and zCompositions. The book was where I learned how to analyze compositional data. I've communicated with several of the dozen-or-so statisticians focused on CoDa, and have a couple of dozen papers, book chapters, and proceedings of the triennial conferences on compositional data analyses. I've also written a monograph that demonstrates how CoDA models applied to benthic macroinvertebrate functional feeding groups can assess water quality and be used to set standards. Thanks, Rich __ 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] Displaying Compositional Data With 46 Parts
Hi Rich, I think this is more a technical question for the subject matter experts than for R-help if I am understanding the question correctly. That said, there seems to be an R package called compositional and a corresponding book http://www.springer.com/us/book/9783642368080 that may help. If nothing else the names of the author(s) of the package or book or the references may give you some leads on some decent papers. It, probably, would help anyone with some expertise in the area to be able to look at some of your data. Have a look at ?dput and perhaps a glance at http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example and http://adv-r.had.co.nz/Reproducibility.html for some suggestions on how to do this. John Kane Kingston ON Canada -Original Message- From: rshep...@appl-ecosys.com Sent: Fri, 17 Jul 2015 06:16:44 -0700 (PDT) To: r-help@r-project.org Subject: [R] Displaying Compositional Data With 46 Parts The compositional data have been divided into two data frames: 46 response variables (the compositional components) and 5 explanatory variables. There are 209 observations of each. With no experience analyzing large compositions with so many parts your advice on how to plot and report results of analyzing these is needed. Matrix plots (scatter, ternary, etc.) of so many parts would be too small when printed on a page to be readable. Haven't found a compositional data water chemistry publication with so many component parts that could be used as a basis for learning how to work with such data sets. Advice and suggestions needed. TIA, Rich __ 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. FREE ONLINE PHOTOSHARING - Share your photos online with your friends and family! Visit http://www.inbox.com/photosharing to find out more! __ 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] Displaying Compositional Data With 46 Parts
The compositional data have been divided into two data frames: 46 response variables (the compositional components) and 5 explanatory variables. There are 209 observations of each. With no experience analyzing large compositions with so many parts your advice on how to plot and report results of analyzing these is needed. Matrix plots (scatter, ternary, etc.) of so many parts would be too small when printed on a page to be readable. Haven't found a compositional data water chemistry publication with so many component parts that could be used as a basis for learning how to work with such data sets. Advice and suggestions needed. TIA, Rich __ 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] TukeyHSD troubles
Bart, I want to thank you for your code. I was having similar problems as Amy, even after setting my numeric variable as a factor using as.factor(). I used is.factor() to confirm and received the answer as TRUE from R; however after running the TukeyHSD() my set factor in my aov() was not read properly. I used your code and it worked perfectly! Thank you! Tracy -- View this message in context: http://r.789695.n4.nabble.com/TukeyHSD-troubles-tp1570205p470.html Sent from the R help mailing list archive at Nabble.com. __ 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] Hausman Test trouble - plm
Hi there. I am a student / very fresh R user who is currently having some issues running the procedure for a Hausman test in R. The head for my data sheet named data looks like this: BirdSeason Gully Grouping Food Habitat.Type 1 83 111 0.152 2 47 111 0.091 3 47 111 1.343 4 47 111 0.091 5 51 113 0.152 6 51 113 0.152 The code to run the test looks like this: library(plm) result=read.csv(H:/data,header=T,sep=,,stringsAsFactors=F) wi=plm(Habitat.Type~Season+Gully+Grouping+Food, data = result, index =c(Bird),model=within) re=plm(Habitat.Type~Season+Gully+Grouping+Food, data = result, index =c(Bird),model=random) phtest(wi, re) The reasoning behind this format is that Habitat.Type is the dependant variable. Season, Gully, Grouping and Food are the independant variables. Bird (ID) is the index as it is the random effect. The error message I am getting is this: duplicate couples (time-id) Error in pdim.default(index[[1]], index[[2]]) From what I have been able to figure out this may be because I have multiple identical observations. The problem is that I do not want to remove these multiple identical observations, as that is a large part of my data. My question to you - am I doing anything wrong? Is there a work around for the duplicate error that I am getting without removing identical observations? Thanks so much for your help. -- View this message in context: http://r.789695.n4.nabble.com/Hausman-Test-trouble-plm-tp4709990.html Sent from the R help mailing list archive at Nabble.com. __ 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] User defined function within a formula
Thanks Bill for your quick reply.
I tried your solution and it did work for the simple user defined function
xploly. But when I try with other function, it gave me error again:
OPoly-function(x,degree=1,weight=1){
weight=round(weight,0)# weight need to be integer
if(length(weight)!=length(x))weight=rep(1,length(x))
p=poly(4*(rep(x,weight)-mean(range(x)))/diff(range(x)),degree)
Z-(t(t(p[cumsum(weight),])*sqrt(attr(p,coefs)$norm2[-seq(2)]))[,degree])
class(Z)-OPoly;Z
}
##this OPoly is an FORTRAN orthogonal polynomial routine, it first maps the
x to range[-2,2] then do QR, then return the results with sqrt(norm2).
Comparing with poly, this transformation will make the model coefficients
within a similar range as other variables, the R poly routine will usually
give you a very large coefficients. I did not find such routine in R, so I
have to define this as user defined function.
###
I also have following function as you suggested:
makepredictcall.OPoly-function(var,call)
{
if (is.call(call)) {
if (identical(call[[1]], quote(OPoly))) {
if (!is.null(tmp - attr(var, coefs))) {
call$coefs - tmp
}
}
}
call
}
But I still got error for following:
g3=glm(lot1 ~ log(u) + OPoly(u,1), data = clotting, family = Gamma)
predict(g3,dc)Error in poly(4 * (rep(x, weight) -
mean(range(x)))/diff(range(x)), degree) :
missing values are not allowed in 'poly'
I thought it might be due to the /diff(range(x) in the function. But even
I remove that part, it will still give me error. Any idea?
Many thanks in advance.
Alex
On Thu, Jul 16, 2015 at 2:09 PM, William Dunlap wdun...@tibco.com wrote:
Read about the 'makepredictcall' generic function. There is a method,
makepredictcall.poly(), for poly() that attaches the polynomial
coefficients
used during the fitting procedure to the call to poly() that predict()
makes.
You ought to supply a similar method for your xpoly(), and xpoly() needs
to return an object of a a new class that will cause that method to be used.
E.g.,
xpoly - function(x,degree=1,...){ ret - poly(x,degree=degree,...);
class(ret) - xpoly ; ret }
makepredictcall.xpoly - function (var, call)
{
if (is.call(call)) {
if (identical(call[[1]], quote(xpoly))) {
if (!is.null(tmp - attr(var, coefs))) {
call$coefs - tmp
}
}
}
call
}
g2 - glm(lot1 ~ log(u) + xpoly(u,1), data = clotting, family = Gamma)
predict(g2,dc)
# 1 2 3 4
5
#-0.01398928608 -0.01398928608 -0.01398928608 -0.01398928608
#-0.01398928608
# 6 7 8 9
#-0.01398928608 -0.01398928608 -0.01398928608 -0.01398928608
You can see the effects of makepredictcall() in the 'terms' component
of glm's output. The 'variables' attribute of it gives the original
function
calls and the 'predvars' attribute gives the calls to be used for
prediction:
attr(g2$terms, variables)
list(lot1, log(u), xpoly(u, 1))
attr(g2$terms, predvars)
list(lot1, log(u), xpoly(u, 1, coefs = list(alpha = 40, norm2 = c(1,
9, 8850
Bill Dunlap
TIBCO Software
wdunlap tibco.com
On Thu, Jul 16, 2015 at 12:35 PM, Kunshan Yin yinkuns...@gmail.com
wrote:
Hello, I have a question about the formula and the user defined function:
I can do following:
###Case 1:
clotting - data.frame(
+ u = c(5,10,15,20,30,40,60,80,100),
+ lot1 = c(118,58,42,35,27,25,21,19,18),
+ lot2 = c(69,35,26,21,18,16,13,12,12))
g1=glm(lot1 ~ log(u) + poly(u,1), data = clotting, family = Gamma)
dc=clotting
dc$u=1
predict(g1,dc)
1 2 3 4 5
6 7 8 9
-0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929 -0.01398929
-0.01398929 -0.01398929 -0.01398929
However, if I just simply wrap the poly as a user defined function ( in
reality I would have my own more complex function) then I will get error:
###Case 2:
xpoly-function(x,degree=1){poly(x,degree)}
g2=glm(lot1 ~ log(u) + xpoly(u,1), data = clotting, family = Gamma)
predict(g2,dc)
Error in poly(x, degree) :
'degree' must be less than number of unique points
It seems that the predict always treat the user defined function in the
formula with I(). My question is how can I get the results for Case2
same
as case1?
Anyone can have any idea about this?
Thank you very much.
Alex
[[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
Re: [R] installation of package ‘devtools’ had non-zero exit status in a POWERPC
I'm stuck on this and i don't know how to install Curl without any problem. Anyone? -- View this message in context: http://r.789695.n4.nabble.com/installation-of-package-devtools-had-non-zero-exit-status-in-a-POWERPC-tp4709687p4709994.html Sent from the R help mailing list archive at Nabble.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] Hausman Test trouble - plm
Might have just solved my own problem team! I assumed that the issue here was the replicated samples, and so added a column and gave a number to each replicate. R seemed to like this and was happy to run the test! A significant result tells me that the fixed effects model is the most preferable model to explain the variation seen in my data. Unless I am doing/assuming something wrong here that you can see then I might well have solved my own problem. Let me know if you have any thoughts :) Cheers -- View this message in context: http://r.789695.n4.nabble.com/Hausman-Test-trouble-plm-tp4709990p4709992.html Sent from the R help mailing list archive at Nabble.com. __ 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] gridPackage not available in r
Dear Responsible, Hello! I am a PhD student at Université Catholique de Louvain, I contact you because I have some difficulties to install neuralnet package in R for my research. However, to install this package, we need two packages MASS and grid. The grid package is not available in my R version. I use R 3.1.1 version for my research. Please help me. I need to develop neural netwok in my research. For example, I have write very well my script but, if I run, the problem is grid package not availbable. install.packages(grid) Installing package into ‘C:/Users/issoufou/Documents/R/win-library/3.2’ (as ‘lib’ is unspecified) Warning in install.packages : package ‘grid’ is not available (for R version 3.2.0) Best regards *Issoufou Ouedraogo* *PhD Student* *Earth and Life Insitute/ Environmental Sciences* *Université Catholique de Louvain**/Belgique* *Croix du sud 2, bte 1, B-1348, Louvain la Neuve,Belgique* *Tel: +32 (0) 10/47.37.19 %2B32%20%280%29%2010%2F47.37.19* *Fax: +32 (0) 10/47.47.45 %2B32%20%280%29%2010%2F47.47.45* [[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] User defined function within a formula
Thank you very much. It worked. I think I need to digest this further
later. Thanks again for the help.
On Thu, Jul 16, 2015 at 4:51 PM, William Dunlap wdun...@tibco.com wrote:
This might do what you want:
OPoly - function(x, degree=1, weight=1, coefs=NULL, rangeX=NULL){
weight - round(weight,0)# weight need to be integer
if(length(weight)!=length(x)) {
weight - rep(1,length(x))
}
if (is.null(rangeX)) {
rangeX - range(x)
}
p - poly(4*(rep(x,weight)-mean(rangeX))/diff(rangeX), degree=degree,
coefs=coefs)
# why t(t(...))? That strips the attributes.
Z - t( t(p[cumsum(weight),]) * sqrt(attr(p,coefs)$norm2[-seq(2)]) )[,
degree, drop=FALSE]
class(Z) - OPoly
attr(Z, coefs) - attr(p, coefs)
attr(Z, rangeX) - rangeX
Z
}
makepredictcall.OPoly-function(var,call)
{
if (is.call(call)) {
if (identical(call[[1]], quote(OPoly))) {
if (!is.null(tmp - attr(var, coefs))) {
call$coefs - tmp
}
if (!is.null(tmp - attr(var, rangeX))) {
call$rangeX - tmp
}
call$weight - 1 # weight not relevant in predictions
}
}
call
}
d - data.frame(Y=1:8, X=log(1:8), Weight=1:8)
fit - lm(data=d, Y ~ OPoly(X, degree=2, weight=Weight))
predict(fit)[c(3,8)]
predict(fit, newdata=data.frame(X=d$X[c(3,8)])) # same result
Bill Dunlap
TIBCO Software
wdunlap tibco.com
On Thu, Jul 16, 2015 at 4:39 PM, William Dunlap wdun...@tibco.com wrote:
OPoly-function(x,degree=1,weight=1){
weight=round(weight,0)# weight need to be integer
if(length(weight)!=length(x))weight=rep(1,length(x))
p=poly(4*(rep(x,weight)-mean(range(x)))/diff(range(x)),degree)
Z-(t(t(p[cumsum(weight),])*sqrt(attr(p,coefs)$norm2[-
seq(2)]))[,degree])
class(Z)-OPoly;Z
}
You need to make OPoly to have optional argument(s) that give
the original-regressor-dependent information to OPoly and then
have it return, as attributes, the value of those arguments.
makepredictcall
will take the attributes and attach them to the call in predvars so
predict uses values derived from the original regressors, not value
derived
from the data to be predicted from.
Take a look at a pair like makepredictcall.scale() and scale() for an
example:
scale has optional arguments 'center' and 'scale' that it returns as
attributes
and makepredictcall.scale adds those to the call to scale that it is
given.
Thus when you predict, the scale and center arguments come from the
original data, not from the data you are predicting from.
Bill Dunlap
TIBCO Software
wdunlap tibco.com
On Thu, Jul 16, 2015 at 3:43 PM, Kunshan Yin yinkuns...@gmail.com
wrote:
Thanks Bill for your quick reply.
I tried your solution and it did work for the simple user defined
function xploly. But when I try with other function, it gave me error again:
OPoly-function(x,degree=1,weight=1){
weight=round(weight,0)# weight need to be integer
if(length(weight)!=length(x))weight=rep(1,length(x))
p=poly(4*(rep(x,weight)-mean(range(x)))/diff(range(x)),degree)
Z-(t(t(p[cumsum(weight),])*sqrt(attr(p,coefs)$norm2[-seq(2)]))[,degree])
class(Z)-OPoly;Z
}
##this OPoly is an FORTRAN orthogonal polynomial routine, it first maps
the x to range[-2,2] then do QR, then return the results with sqrt(norm2).
Comparing with poly, this transformation will make the model coefficients
within a similar range as other variables, the R poly routine will usually
give you a very large coefficients. I did not find such routine in R, so I
have to define this as user defined function.
###
I also have following function as you suggested:
makepredictcall.OPoly-function(var,call)
{
if (is.call(call)) {
if (identical(call[[1]], quote(OPoly))) {
if (!is.null(tmp - attr(var, coefs))) {
call$coefs - tmp
}
}
}
call
}
But I still got error for following:
g3=glm(lot1 ~ log(u) + OPoly(u,1), data = clotting, family = Gamma)
predict(g3,dc)Error in poly(4 * (rep(x, weight) -
mean(range(x)))/diff(range(x)), degree) :
missing values are not allowed in 'poly'
I thought it might be due to the /diff(range(x) in the function. But
even I remove that part, it will still give me error. Any idea?
Many thanks in advance.
Alex
On Thu, Jul 16, 2015 at 2:09 PM, William Dunlap wdun...@tibco.com
wrote:
Read about the 'makepredictcall' generic function. There is a method,
makepredictcall.poly(), for poly() that attaches the polynomial
coefficients
used during the fitting procedure to the call to poly() that predict()
makes.
You ought to supply a similar method for your xpoly(), and xpoly()
needs to return an object of a a new class that will cause that method to
be used.
E.g.,
xpoly - function(x,degree=1,...){ ret - poly(x,degree=degree,...);
class(ret) - xpoly ; ret }
makepredictcall.xpoly - function (var, call)
{
if (is.call(call)) {
if (identical(call[[1]],
Re: [R] installation of package ‘devtools’ had non-zero exit status in a POWERPC
That is an operating-system configuration question, not a question about R. There are many OSs out there... please find a forum with users of your OS in which to pursue this question. --- Jeff NewmillerThe . . Go Live... DCN:jdnew...@dcn.davis.ca.usBasics: ##.#. ##.#. Live Go... Live: OO#.. Dead: OO#.. Playing Research Engineer (Solar/BatteriesO.O#. #.O#. with /Software/Embedded Controllers) .OO#. .OO#. rocks...1k --- Sent from my phone. Please excuse my brevity. On July 17, 2015 4:04:14 AM PDT, aidaph ap...@alumnos.unican.es wrote: I'm stuck on this and i don't know how to install Curl without any problem. Anyone? -- View this message in context: http://r.789695.n4.nabble.com/installation-of-package-devtools-had-non-zero-exit-status-in-a-POWERPC-tp4709687p4709994.html Sent from the R help mailing list archive at Nabble.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. __ 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] Displaying Compositional Data With 46 Parts
See in-line John Kane Kingston ON Canada -Original Message- From: rshep...@appl-ecosys.com Sent: Fri, 17 Jul 2015 07:33:43 -0700 (PDT) To: r-help@r-project.org Subject: Re: [R] Displaying Compositional Data With 46 Parts On Fri, 17 Jul 2015, John Kane wrote: I think this is more a technical question for the subject matter experts than for R-help if I am understanding the question correctly. John, I agree completely. Unfortunately, there is no R SIG devoted to CoDA, nor any other mail list or Web forum that I've been able to find. There is a CoDA Web site but no active forum. That said, there seems to be an R package called compositional and a corresponding book http://www.springer.com/us/book/9783642368080 that may help. If nothing else the names of the author(s) of the package or book or the references may give you some leads on some decent papers. In addition to compositions, there's robCompositions and zCompositions. The book was where I learned how to analyze compositional data. I've communicated with several of the dozen-or-so statisticians focused on CoDa, and have a couple of dozen papers, book chapters, and proceedings of the triennial conferences on compositional data analyses. I've also written a monograph that demonstrates how CoDA models applied to benthic macroinvertebrate functional feeding groups can assess water quality and be used to set standards. Then it sounds like you are one of the experts. Do whatever you think appropriate and either set the standard for future research or get enough feedback to do even better next time. :) Thanks, Rich __ 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. Send your photos by email in seconds... TRY FREE IM TOOLPACK at http://www.imtoolpack.com/default.aspx?rc=if3 Works in all emails, instant messengers, blogs, forums and social networks. __ 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] powerTransform warning message?
Thank you so much for the explanation. That was very helpful! :-) Thanks! Brittany On Jul 16, 2015, at 6:16 PM, John Fox j...@mcmaster.ca wrote: Dear Brittany, On Thu, 16 Jul 2015 17:35:38 -0600 Brittany Demmitt demmi...@gmail.com wrote: Hello, I have a series of 40 variables that I am trying to transform via the boxcox method using the powerTransfrom function in R. I have no zero values in any of my variables. When I run the powerTransform function on the full data set I get the following warning. Warning message: In sqrt(diag(solve(res$hessian))) : NaNs produced However, when I analyze the variables in groups, rather than all 40 at a time I do not get this warning message. Why would this be? And does this mean this warning is safe to ignore? No, it is not safe to ignore the warning, and the problem has nothing to do with non-positive values in the data -- when you say that there are no 0s in the data, I assume that you mean that the data values are all positive. The square-roots of the diagonal entries of the Hessian at the (pseudo-) ML estimates are the SEs of the estimated transformation parameters. If the Hessian can't be inverted, that usually implies that the maximum of the (pseudo-) likelihood isn't well defined. This isn't surprising when you're trying to transform as many as 40 variables at a time to multivariate normality. It's my general experience that people often throw their data into the Box-Cox black box and hope for the best without first examining the data, and, e.g., insuring a reasonable ratio of maximum/minimum values for each variable, checking for extreme outliers, etc. Of course, I don't know that you did that, and it's perfectly possible that you were careful. I would like to add that all of my lambda values are in the -5 to 5 range. I also get different lambda values when I analyze the variables together versus in groups. Is this to be expected? Yes. It's very unlikely that both are right. If, e.g., the variables are multivariate normal within groups then their marginal distribution is a mixture of multivariate normals, which almost surely isn't itself normal. I hope this helps, John John Fox, Professor McMaster University Hamilton, Ontario, Canada http://socserv.mcmaster.ca/jfox/ Thank you so much! Brittany [[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.
Re: [R] removing the columns with 0 or NA or 1or NA or 2 or NA
Hi Lida, I think that your matrix is actually a data frame, so try this: mat[,sapply(mat,function(x) var(x,na.rm=TRUE)0)] Jim On Fri, Jul 17, 2015 at 12:58 AM, Lida Zeighami lid.z...@gmail.com wrote: I have ma matrix which its elements are NA,0,1,2 ! I got my answer bout removing the columns with 0 or NA or both values but now I want to add additional condition for deleting the columns! I have to delete the columns which contain the same value. delete the columns with NA or 0 or both and the columns with NA or 1 or both and the column with NA or 2 or both (I should keep the columns which have variation in their values)! I use this code but didn't work properly: mat_nonNA- mat[, !apply((is.na(mat) | mat == 0) (is.na(mat) | mat==1) ( is.na(mat) | mat==2), 2, all)] mat 1:1105901701:110888172 1:110906406 1:110993854 1:110996710 1:44756 A05363 1 1 1 2 NA 0 A05370 0 1 0NA 0 NA A05380 1 NA 2NA NA 0 A05397 01 0NA 0 2 A05400 21 0 2 0 0 A05426 0 NA NA NA 0 1 my out put should be like below: 1:110590170 1:110906406 1:44756 A05363 1 1 0 A05370 0 0 NA A05380 1 2 0 A05397 0 0 2 A05400 2 0 0 A05426 0 NA 1 Thanks for your help [[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] Release of R 3.2.2 scheduled for August 14
We intend to have a patch release on August 14, nickname will be Fire Safety. The detailed schedule will be made available via developer.r-project.org as usual. -- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Office: A 4.23 Email: pd@cbs.dk Priv: pda...@gmail.com ___ r-annou...@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-announce __ 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] Displaying Compositional Data With 46 Parts
On Fri, 17 Jul 2015, Bert Gunter wrote: I believe John Aitchison's book and papers are the authoritative basic resources. Have you read them? Bert, Yes, I have. The problem is that the support of the distributions are (hyper)simplexes, not Euclidean space, due to the requirement that the proportions must sum to 1. This means that complicated animals like Dirichelet distributions must be used to model populations, and the sampling theory is therefore specialized. It's difficult for most folks to get their heads around this. That's true, When I read the math I move my lips and follow along with a finger. :-) My question is focused on presentation of graphic presentation of the data, such as a matrix of ternary diagrams that show the distribution of the response variables to the explanatory variables. The analysis of benthic macroinvertebrate functional feeding groups has 5 response variables which resulted in a 5-by-5 ternary diagram matrix. Anything larger than that would require the eyesight of a teenager to see any details. I'll keep searching the literature for a suitable example. Perhaps a CoDA SIG will develop on within the R mail list ecosystem in the not-too-distant future. Thanks, Rich __ 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] Displaying Compositional Data With 46 Parts
On Fri, 17 Jul 2015, John Kane wrote: Then it sounds like you are one of the experts. Do whatever you think appropriate and either set the standard for future research or get enough feedback to do even better next time. :) John, Far from an expert, but becoming more capable with each project. Someone, somewhere, has dealt successfully with this issue in the past. So I need to keep inquiring and hoping that e-mail addresses are still valid and that I get responses. :-) Compositional data are the rule for public (e.g., economic, political, societal) and environmental chemical data, including some biological data. It seems to be unknown outside of those in academia focused on developing the mathematical theory and applying the models to their research. Regards, Rich __ 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] gridPackage not available in r
On Jul 17, 2015, at 7:21 AM, Issoufou Ouedraogo wrote: Dear Responsible, Hello! I am a PhD student at Université Catholique de Louvain, I contact you because I have some difficulties to install neuralnet package in R for my research. However, to install this package, we need two packages MASS and grid. The grid package is not available in my R version. I use R 3.1.1 version for my research. Please help me. I need to develop neural netwok in my research. For example, I have write very well my script but, if I run, the problem is grid package not availbable. install.packages(grid) The grid package is part of the standard R installation. It is not available on CRAN as a separate package. Just do this: library(grid) # you should get a message help(pac=grid) -- David. Installing package into ‘C:/Users/issoufou/Documents/R/win-library/3.2’ (as ‘lib’ is unspecified) Warning in install.packages : package ‘grid’ is not available (for R version 3.2.0) Best regards *Issoufou Ouedraogo* *PhD Student* *Earth and Life Insitute/ Environmental Sciences* *Université Catholique de Louvain**/Belgique* *Croix du sud 2, bte 1, B-1348, Louvain la Neuve,Belgique* *Tel: +32 (0) 10/47.37.19 %2B32%20%280%29%2010%2F47.37.19* *Fax: +32 (0) 10/47.47.45 %2B32%20%280%29%2010%2F47.47.45* [[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. David Winsemius Alameda, CA, USA __ 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] OPTIMX: non-finite finite-difference value [26] and scaling problem
Hello,I am running a nonlinear GMM using the optimx wrapper. I am trying to estimate 37 variables however and my code for the optimx is: nlgmm = optimx(par=b0, fn=obj,method = BFGS, itnmax=1, control=list(follow.on = TRUE,kkt=FALSE,starttests=TRUE,save.failures=TRUE, trace=0)) My staring values are from Ordinary least squares estimates (OLS) and they are: b0 - c(-2.00658,-0.04373,0.19079,0.34652,-0.36814,0.21284,-0.24369,0.64622,0.22927,0.29431,0.19547,0.80614,18.8398,0.5928,3.1375,0.4301,-0.4937,2.2016,31.5203,0.6171,1.0206,1.7830,-0.4421,11.0076,-0.03305,0.17087,0.44794,0.17488,0.10781,-0.50747,-0.04563,0.17030,0.41792,0.17526,0.99734,-17.2996,-41.9359) I got the following error: Error in optim(par = par, fn = ufn, gr = ugr, lower = lower, upper = upper, : non-finite finite-difference value [26] I also got an error about scaling which is as follows: Parameters or bounds appear to have differentscalings.This can cause poor performance in optimization. Itis important for derivative free methods like BOBYQA, UOBYQA, NEWUOA. any help will be greatly appreciated. Best Regards [[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] Nagelkerke Pseudo R-squared
Dear R-Experts, I have fitted an ordinal logistic regression with just 1 explanatory variable for the reproducible example here below. Everything is working, now I try to calculate the Nagelkerke Pseudo R-squared. I have found a package BaylorEdPsych providing many Pseudo R-squared, but the example shown in the package is for GLM (binary logistic regression) not for ordinal logistic regression. How can I calculate the Nagelkerke Pseudo R-squared for ordinal logistic regression ? Many thanks as usual for your precious help. Reproducible example : install.packages(MASS) library(MASS) a=factor(c(tres grand, grand, petit,petit,tres grand,grand,petit,petit,tres grand,grand)) b=c(homme, homme, femme, femme, femme, homme, homme, homme, femme, femme) m - polr(a ~ b, Hess=TRUE) summary(m) __ 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] Warning message with maxLik()
Dear All,
I'm trying to get the MLe for a certain distribution using maxLik ()
function. I wrote the log-likelihood function as follows:
theta -vector(mode = numeric, length = 3)
r- 17
n -30
T-c(7.048,0.743,2.404,1.374,2.233,1.52,23.531,5.182,4.502,1.362,1.15,1.86,1.692,11.659,1.631,2.212,5.451)
C-
c(0.562,5.69,12.603,3.999,6.156,4.004,5.248,4.878,7.122,17.069,23.996,1.538,7.792)
# The loglik. func.
loglik - function(param) {
theta[1]- param[1]
theta[2]- param[2]
theta[3]- param[3]
l-(r*log(theta[3]))+(r*log(theta[1]+theta[2]))+(n*theta[3]*log(theta[1]))+(n*theta[3]*log(theta[2]))+
(-1*(theta[3]+1))*sum(log((T*(theta[1]+theta[2]))+(theta[1]*theta[2])))+
(-1*theta[3]*sum(log((C*(theta[1]+theta[2]))+(theta[1]*theta[2]
return(l)
}
then, I evaluated it at theta- c(40,50,2)
v-loglik(param=theta)
v
[1] -56.66653
I used this same log-likelihood function, once with analytic gradient and
another time with numerical one, with the maxLik function, and in both
cases I got the same 50 warning messages and an MLE which is completely
unrealistic as per my applied example.
a - maxLik(loglik, gradlik, hesslik, start=c(40,50,2))
where gradlik and hesslik are the analytic gradient and Hessian matrix,
respectively, given by:
U - vector(mode=numeric,length=3)
gradlik-function(param = theta,n, T,C)
{
U - vector(mode=numeric,length=3)
theta[1] - param[1]
theta[2] - param[2]
theta[3] - param[3]
r- 17
n -30
T-c(7.048,0.743,2.404,1.374,2.233,1.52,23.531,5.182,4.502,1.362,1.15,1.86,1.692,11.659,1.631,2.212,5.451)
C-
c(0.562,5.69,12.603,3.999,6.156,4.004,5.248,4.878,7.122,17.069,23.996,1.538,7.792)
U[1]- (r/(theta[1]+theta[2]))+((n*theta[3])/theta[1])+(
-1*(theta[3]+1))*sum((T+theta[2])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+
(-1*(theta[3]))*sum((C+theta[2])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
U[2]-(r/(theta[1]+theta[2]))+((n*theta[3])/theta[2])+
(-1*(theta[3]+1))*sum((T+theta[1])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+
(-1*(theta[3]))*sum((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
U[3]-(r/theta[3])+(n*log(theta[1]*theta[2]))+
(-1)*sum(log((T*(theta[1]+theta[2]))+(theta[1]*theta[2])))+(-1)*sum(log((C*(theta[1]+theta[2]))+(theta[1]*theta[2])))
return(U)
}
hesslik-function(param=theta,n,T,C)
{
theta[1] - param[1]
theta[2] - param[2]
theta[3] - param[3]
r- 17
n -30
T-c(7.048,0.743,2.404,1.374,2.233,1.52,23.531,5.182,4.502,1.362,1.15,1.86,1.692,11.659,1.631,2.212,5.451)
C-
c(0.562,5.69,12.603,3.999,6.156,4.004,5.248,4.878,7.122,17.069,23.996,1.538,7.792)
G- matrix(nrow=3,ncol=3)
G[1,1]-((-1*r)/((theta[1]+theta[2])^2))+((-1*n*theta[3])/(theta[1])^2)+
(theta[3]+1)*sum(((T+theta[2])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))^2)+(
theta[3])*sum(((C+theta[2])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))^2)
G[1,2]-((-1*r)/((theta[1]+theta[2])^2))+
(theta[3]+1)*sum(((T)/((theta[1]+theta[2])*T+(theta[1]*theta[2])))^2)+
(theta[3])*sum(((C)/((theta[1]+theta[2])*C+(theta[1]*theta[2])))^2)
G[2,1]-G[1,2]
G[1,3]-(n/theta[1])+(-1)*sum(
(T+theta[2])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+(-1)*sum((C+theta[2])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
G[3,1]-G[1,3]
G[2,2]-((-1*r)/((theta[1]+theta[2])^2))+((-1*n*theta[3])/(theta[2])^2)+
(theta[3]+1)*sum(((T+theta[1])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))^2)+(
theta[3])*sum(((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))^2)
G[2,3]-(n/theta[2])+(-1)*sum((T+theta[1])/((theta[1]+theta[2])*T+(theta[1]*theta[2])))+(-1)*sum((C+theta[1])/((theta[1]+theta[2])*C+(theta[1]*theta[2])))
G[3,2]-G[2,3]
G[3,3]-((-1*r)/(theta[3])^2)
return(G)
}
and using numeric gradient and hessian matrix:
a - maxLik(loglik, start=c(40,50,2))
Warning messages:
1: In log(theta[3]) : NaNs produced
2: In log(theta[1] + theta[2]) : NaNs produced
3: In log(theta[1]) : NaNs produced
4: In log((T * (theta[1] + theta[2])) + (theta[1] * theta[2])) : NaNs
produced
5: In log((C * (theta[1] + theta[2])) + (theta[1] * theta[2])) : NaNs
produced
6: In log(theta[3]) : NaNs produced
7: In log(theta[1] + theta[2]) : NaNs produced
and so on…..
I don't know why I get these 50 warnings although:
1- The inputs of the log() function are strictly positive.
2- When I evaluated the log-likelihood fuction at the very begining it gave
me a number(which is -56.66) and not (NAN).
I've also tried to:
1- Reparamtrize my model using lamda(i)= log(theta(i)), for i=1,2,3, so
that it may solve the problem, but it didn't.
2- I've used the comparederivitive() function, and the analytic and numeric
gradients were so close.
Any help please?
Maram Salem
[[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] matching strings in a list
Thank you all very much! -- View this message in context: http://r.789695.n4.nabble.com/matching-strings-in-a-list-tp4709967p4710015.html Sent from the R help mailing list archive at Nabble.com. __ 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] modifying a package installed via GitHub
Hi Folks, I am working with a package installed via GitHub that I would like to modify. However, I am not sure how I would go about loading a 'local' version of the package after I have modified it, and whether that process would including uninstalling the original unmodified package (and, conversely, how to uninstall my local, modified version if I wanted to go back to the unmodified version available on GitHub). Any advice would be appreciated. Thanks, Steve -- View this message in context: http://r.789695.n4.nabble.com/modifying-a-package-installed-via-GitHub-tp4710016.html Sent from the R help mailing list archive at Nabble.com. __ 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] Displaying Compositional Data With 46 Parts
Hi Rich, Being in a position of relative ignorance on this topic, I'll offer some suggestions that may well be useless. You mention ternary diagrams, which use position to represent compositional proportions. These will not scale up to 46 values in any way that I can imagine. If you want to display relative concentration or the like, differentiate the components and display numeric information about each component, you may be looking for a variant of the Hinton diagram. This is a bit like a heatmap where the squares are of different sizes, representing relative numeric values. In the Hinton diagram, the colors represent sign (+-), but you would probably want more complex coloring. Finally, identifying labels and/or values could be displayed on each cell of the matrix. It would not be too difficult to program something like this, so if this idea is not completely useless, let me know. It is of course possible to go the interactive route and produce a play with me display if necessary. 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.
Re: [R] Nagelkerke Pseudo R-squared
On Jul 17, 2015, at 4:33 PM, varin sacha wrote: Dear R-Experts, I have fitted an ordinal logistic regression with just 1 explanatory variable for the reproducible example here below. Everything is working, now I try to calculate the Nagelkerke Pseudo R-squared. I have found a package BaylorEdPsych providing many Pseudo R-squared, but the example shown in the package is for GLM (binary logistic regression) not for ordinal logistic regression. How can I calculate the Nagelkerke Pseudo R-squared for ordinal logistic regression ? polr-objects have a deviance node. If this has statistical value (which I have some doubts regarding) then just apply the usual formula to compare nested models. -- David. Many thanks as usual for your precious help. Reproducible example : install.packages(MASS) library(MASS) a=factor(c(tres grand, grand, petit,petit,tres grand,grand,petit,petit,tres grand,grand)) b=c(homme, homme, femme, femme, femme, homme, homme, homme, femme, femme) m - polr(a ~ b, Hess=TRUE) summary(m) __ 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. David Winsemius Alameda, CA, USA __ 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.
