title: Confidence intervals for a polynomial

text: I have a random variable $Z$ which takes values in the nonnegative integers $\\{ 0,1,2,\dots \\}$, call the probabilities for each outcome $z_k:=P[Z=k]$. I can sample from $Z$'s distribution independently and cheaply; I currently have a sample size of $2^{28}$. It appears that $z_0\approx 0.24, z_1\approx 0.18,\dots$, with roughly exponential decay. I have a sequence of quadratic forms with positive coefficients:   * $Q_0(z_0) = \frac14 z_0^2$   * $Q_1(z_0,z_1) = \frac 12 {z_0 z_1}$   * ...   * $Q_7(z_0,z_1,\dots,z_7) = \frac{1}{8} \left(2 z_0 z_1+3 z_2 z_1+4 z_4 z_1+4 z_6 z_1+3 z_0 z_3 + \right.$ $\qquad\qquad\qquad\qquad\qquad\qquad \left. +4 z_2 z_3+4 z_3 z_4+4 z_0 z_5+4 z_2 z_5+4 z_0 z_7\right)$   * ... What I would like to have is a confidence interval for the $Q_i$'s that is less than $10^{-4}$ wide, but I'll take whatever I can get. I have rigorous bounds on the $z_i$, and since the coefficients of the $Q$'s are all positive, it is straightforward to turn these into rigorous bounds for the $Q$'s. But I don't know how to do this correctly with confidence intervals. * * * What's this about? I found a bizarre phenomenon in number theory, and I know how to prove that it really happens, but actually doing so will require some programming effort on my part and a considerable amount of time on our local cluster. Before I invest that time and clog up our machine, I'd like to be more certain than I am that the phenomenon is real. I want to quantify the reasonableness of my claim that $Q_7<Q_6$ and $Q_7<Q_8$. My estimates indicate that $Q_6-Q_7$ is around $5\cdot 10^{-4}$, which is why I wanted CIs at that resolution. Fix a large integer $n$, and let $A$ be a uniformly chosen subset of $\\{1,2,\dots,n\\}$ (that is, each particular subset has probability $2^{-n}$ of being chosen). Let $Q_k(n)$ be the probability that exactly $k$ of the numbers from $\\{2,3,\dots,2n\\}$ cannot be written as a sum of two elements of $A$; let $Q_k = \lim_n Q_k(n)$. It's a little tricky to prove, but those limits exist and $\sum_{k} Q_k =1$. Now it's no surprise that $Q_0$ is small, and as $k$ increases $Q_k$ increases, has a peak and then decays exponentially. The bizarre part is that there is a bias against 7. That is, experimentally $Q_7< Q_6$ and $Q_7<Q_8$. That is, what wasn't a surprise actually isn't true: the distribution is bimodal. I can express the $Q_i$'s (using some theory) as above without the limit in terms of this other distribution, defined by the $z_i$'s. That's handy because I have a way to rigorously bound the $z_i$'s using, as I mentioned above, some large computations. Also, I have a very large data set for the $Z$ variable.


title: Comparing standard deviation across multiple time horizons

text: I am working on comparing a number of investments and am attempting to do a relative risk comparison but am having troubles ensuring accuracy with the standard deviations. A number of the observations have a shorter time horizon than I would like to compare (some have only a 3,4, or 5 year time horizon). I would like to use a 10-year time horizon so 10 observations in every standard deviation calculation. I believe I can assume a normal distribution for this exercise. Based on this problem, is there a solution and if so, do I need to adjust the shorter standard deviations to ensure accuracy using a longer timeline? This is my first time here and from the looks of it there are some very bright people. I appreciate any help you can provide.


title: Using AIC to select an upper bin for a counting variable

text: I have some patient data and I'm using a logistic model to explore factors that might affect a patient considered severe/not severe for a disease (the DV). Many of the variables are counts; counts of prescriptions in one year, counts of diagnoses codes, GP visits, etc. To take one example, GP visits. Commonly, this is modeled as 0, 1 or 2+ visits a year (so, two dummy variables, each compared with 0 visits). It is never modeled as a continuous variable as risk to outcome is not thought linear, the top of the distribution is often erratic - and physicians hate continuous variables ... ! Using 0/1/2+ feels too simplistic however and I was considering a different method for determining the 'upper bin'. An approach I'm trying is modeling the IV alone against the DV using many bins (e.g. 0/1+, 0/1/2+, 0/1/2/3+) and saving the AICc score for each one. I then plot the score for each model and select the one with the lowest score (with an eye on it being sensible too). Is this a valid approach? Or might there be another way for an inferential model? Thanks, Richard


title: How to handle underdispersion in GLMM (binomial outcome variable)

text: I'm working on the following model in R:                 Generalized linear mixed model fit by maximum likelihood ['glmerMod']       Family: binomial (logit)       Formula: Tooluse ~ Sex + Age + Frequency + Tool.related.skill + (1|Trial) + (1 +                               Frequency|Subjectnumber) + (1 + Tool.related.skill|Frequency/Task)        Data: g4            with Tooluse (yes, no) age (continuous) tool.related.skill (ordinal) trial (1-4) frequency (low, high) task (1-12, nested within frequency. 6 tasks belong to the low frequency group, 6 tasks to high frequency) My research question looks at the effect of the frequency variable on tool use. Testing the model assumptions, I get this output for the test of overdispersion: > overdisp.test (B1NF.FULL) > >> >>        chisq     df    P   dispersion.parameter >>     1 36.68702  141    1      0.2601916 >>   How can I deal with the problem of **underdispersion**? So far I got 3 suggestions (2 of them from one of the authors of the lme4 package): 1) using **mixture/hurdle** models 2) allowing a **negative correlation structure** within groups (which can't be done with lme4 and is harder for GLMMs in general) 3) standard **'quasi-likelihood' approach** , i.e. taking the estimated level of underdispersion and shrinking all the confidence intervals accordingly as a first approach. However, I got warned that the thing to be careful about there is that it has yet to be figured out how quasi-likelihood estimates of 'residual' variance interact with the estimates of the random effects variances I would greatly appreciate any opinions and especially any help on how to implement any of these strategies in R. I feel kind of lost here. Thank you very much in advance for your help.


title: Compute and graph the LDA decision boundary

text: I saw an LDA (linear discriminant analysis) plot with decision boundaries from _The Elements of Statistical Learning_: ![enter image description here](http://i.stack.imgur.com/FBuR7.png) I understand that data are projected onto a lower-dimensional subspace. However, I would like to know how we get the decision boundaries in the original dimension such that I can project the decision boundaries onto a lower-dimensional subspace (likes the black lines in the image above). Is there a formula that I can use to compute the decision boundaries in the original (higher) dimension? If yes, then what inputs does this formula need?


title: R: What do I see in partial dependence plots of gbm and RandomForest?

text: Actually, I thought I had understood what one can show a with partial dependence plot, but using a very simple hypothetical example, I got rather puzzled. In the following chunk of code I generate three independent variables ( _a_ , _b_ , _c_ ) and one dependent variable ( _y_ ) with _c_ showing a close linear relationship with _y_ , while _a_ and _b_ are uncorrelated with _y_. I make a regression analysis with a boosted regression tree using the R package `gbm`:               a <- runif(100, 1, 100)     b <- runif(100, 1, 100)     c <- 1:100 + rnorm(100, mean = 0, sd = 5)     y <- 1:100 + rnorm(100, mean = 0, sd = 5)     par(mfrow = c(2,2))     plot(y ~ a); plot(y ~ b); plot(y ~ c)     Data <- data.frame(matrix(c(y, a, b, c), ncol = 4))     names(Data) <- c("y", "a", "b", "c")     library(gbm)     gbm.gaus <- gbm(y ~ a + b + c, data = Data, distribution = "gaussian")     par(mfrow = c(2,2))     plot(gbm.gaus, i.var = 1)     plot(gbm.gaus, i.var = 2)     plot(gbm.gaus, i.var = 3)      Not surprisingly, for variables _a_ and _b_ the partial dependence plots yield horizontal lines around the mean of _a_. What me puzzles is the plot for variable _c_. I get horizontal lines for the ranges _c_ < 40 and _c_ > 60 and the y-axis is restricted to values close to the mean of _y_. Since _a_ and _b_ are completely unrelated to _y_ (and thus there variable importance in the model is 0), I expected that _c_ would show partial dependence along its entire range instead of that sigmoid shape for a very restricted range of its values. I tried to find information in Friedman (2001) "Greedy function approximation: a gradient boosting machine" and in Hastie et al. (2011) "Elements of Statistical Learning", but my mathematical skills are too low to understand all the equations and formulae therein. Thus my question: What determines the shape of the partial dependence plot for variable _c_? (Please explain in words comprehensible to a non-mathematician!) ADDED on 17th April 2014: While waiting for a response, I used the same example data for an analysis with R-package `randomForest`. The partial dependence plots of randomForest resemble much more to what I expected from the gbm plots: the partial dependence of explanatory variables _a_ and _b_ vary randomly and closely around 50, while explanatory variable _c_ shows partial dependence over its entire range (and over almost the entire range of _y_ ). What could be the reasons for these different shapes of the partial dependence plots in `gbm` and `randomForest`? ![partial plots of gbm and randomForest](http://i.stack.imgur.com/PrlC1.jpg) Here the modified code that compares the plots:               a <- runif(100, 1, 100)     b <- runif(100, 1, 100)     c <- 1:100 + rnorm(100, mean = 0, sd = 5)     y <- 1:100 + rnorm(100, mean = 0, sd = 5)     par(mfrow = c(2,2))     plot(y ~ a); plot(y ~ b); plot(y ~ c)     Data <- data.frame(matrix(c(y, a, b, c), ncol = 4))     names(Data) <- c("y", "a", "b", "c")          library(gbm)     gbm.gaus <- gbm(y ~ a + b + c, data = Data, distribution = "gaussian")          library(randomForest)     rf.model <- randomForest(y ~ a + b + c, data = Data)          x11(height = 8, width = 5)     par(mfrow = c(3,2))     par(oma = c(1,1,4,1))     plot(gbm.gaus, i.var = 1)     partialPlot(rf.model, Data[,2:4], x.var = "a")     plot(gbm.gaus, i.var = 2)     partialPlot(rf.model, Data[,2:4], x.var = "b")     plot(gbm.gaus, i.var = 3)     partialPlot(rf.model, Data[,2:4], x.var = "c")     title(main = "Boosted regression tree", outer = TRUE, adj = 0.15)     title(main = "Random forest", outer = TRUE, adj = 0.85)


title: Can I use a regression model with ANOVA significance greater than 0.05?

text: I have a multiple regression using SPSS. The significance of my model in the ANOVA table is p=0.174 which is >0.05. what does this mean for my model? Can I still use it and proceed in the interpretation of the results or not?


title: Forest plot for meta-analysis displaying the mean ES with and without outliers

text: I am conducting a meta-analysis and am trying to produce a forest plot displaying the mean weighted effect size with and without the outliers. I looked on so many websites and tried a lot of syntaxes, however, didn't really find anything about what I am looking for. Assume there are 10 studies in the meta-analysis. Outlier analyses showed that Study 3 is an outlier leading to an overestimation of the mean weighted ES. Accordingly, I calculate a mean ES with ( _k_ = 10) and without the outlier ( _k_ = 9). I would like to display the studies sorted according to their ES and report both Mean ESs in the forest plot. I would appreciate any help and/or advice. Study3 Study9 Study5 ... Study7 Mean ES ( _k_ = 10) Mean ES ( _k_ = 9) * * * Thank you again for your response. I have been able to produce the graphs I wanted, there were only two issues, and I was wondering if there is an easier way to do it than I did. APA style requires to write statistical terms italic, e.g., _k_ = 10, and at sometimes also subscripted, e.g., u in Hedges _g[u]_. My apologies, the subscript command here seems not to work. I tried this syntax to finish the forest plot but without success as the program would draw the ** and [], as well as the command "italic()" in the graphs. Further, the command digits=2 was successful for the single ES and CIs, however not for the tables left hand side where I have mean age and the mean scores for each of the experimental groups (Note: the DV is continuous). I created a syntax as follows, however, as I said, I was wondering if there is an easier was to do this. The meta-analyses for each of the subgroups have been already calculated.               ### decrease margins so the full space is used     par(mar=c(4,4,1,2))          ### draw a plot     par(font=1)     with(SMetaHits, forest(rma.allHits, xlim=c(-12,  6), ylim=c(-1, 41),                            xlab="Hedges gu", slab=paste(author),                             ilab=cbind(SMetaHits$N, SMetaHits$MYou, SMetaHits$MOld),                             ilab.xpos=c(-6, -4.5, -3.5), ilab.pos=c(2, 4, 4),                            order=order(SMetaHits$FacAgeNum), rows=c(36:30, 24:15, 9:4)))     text(-12, 40, "Author(s), Year [Exp.]", pos=4)          ### add text for the subgroups     par(font=3)     text(-12.1, c(37.5,25.5,10.5), pos=4, c("Young", "Mixed-age", "Old"))      ### This is where it gets complicated and I was hoping that there is an easier way than this               ### add column headings to the plot and draw Y and O subscripted and italic     par(font=3)     text(c(-6), 40, c("N"), pos=2)     text(c(-4.5, -3.5), 40, c("M", "M"), pos=4)     text(c(-4.3, -3.3), 39.7, c("Y", "O"), pos=4, cex=.75) ##Y and O subscripted          ###Command to draw 0 to the mean scores; otherwise e.g., .6 instead of .60.      par(font=1)     text(-3.75, c(36,16,4), "0", pos=2) ###add 0 to MY      text(-2.75, c(7,4), "0", pos=2) ###add 0 to MO      At this stage, there were 2 problems: 1. I chose pos=2 to have all values left handed, which is a good solution only if the values are all in the same direction. As soon as there is one value negative, the values are not exactly underneath each other anymore. There are a lot of trial and errors to get the 0 at the right place, a similar problem I had below with _gu_               ###Command to write Hedges *gu* as required, i.e., gu italic and u subscript.      text(6, 40, "Hedges     [95% CI]", pos=2)       par(font=3)       text(4.93, 40, "g", pos=2)     text(5.02, 39.7, "u", pos=2, cex=.75)      The last command, I noticed only now that _k_ is not italic yet               ### add summary polygons for the three subgroups      par(font=2)     addpoly(res.YoF, row=28.5, cex=1, mlab="FE Model (k = 7)")     addpoly(res.YoFNoOutl, row=27.3, cex=1, mlab="FE Model (k = 6)")     addpoly(res.MxF, row=13.5, cex=1, mlab="FE Model (k = 10)")     addpoly(res.MxFNoOutl, row=12.3, cex=1, mlab="FE Model (k = 9)")     addpoly(res.OlF, row=2.5, cex=1, mlab="FE Model (k = 6)")     addpoly(res.OlFNoOutl, row=1.3, cex=1, mlab="FE Model (k = 5)")}'


title: Is confidence transitive in asociation rules?

text: Given a set of rule such as:   $A \rightarrow B$; $B \rightarrow C$; that satisfy minimum confidence in the context of apriori algorithm meaning: $$\text{Conf}(A \rightarrow B) \geq \text{min. confidence}$$ and $$\text{Conf}(B \rightarrow C) \geq \text{min. confidence}$$ Is it possible to state that $\text{Conf}(A \rightarrow C) \geq \text{min. confidence}$? I believe it's not possible to assert the later given that the relation between $\text{Support}(A \rightarrow B)$ and $\text{Support}(B \rightarrow C)$ doesn't give any information about the $\text{Support}(A \rightarrow C)$ Does this make sense? Note: it's just a book exercise I'm trying to solve(not homework!) Here is the chapter link, it's in page 405 exercise 3.


title: Rice distribution: estimate $\nu$ from 2D-data

text: How can I get a good estimate of parameter $\nu$ of the Rice distribution based on a set of $(x,y)$-coordinates? **Edit** : Given whuber's excellent comment, I'm not looking for an _unbiased_ estimate. Instead, I'd like to know whether there is an estimate that takes advantage of having available the 2D points themselves (not just their Euclidean norm) as well as an estimate of $\sigma$. Background: Given an uncorrelated bivariate normal variable $X$ with mean $\mu \neq 0$ with covariance matrix $\Sigma = \sigma^{2}I$, the Euclidean norm $R = \sqrt{X'X}$ (the radius, or distance to the origin) follows a Rice distribution with parameters $\nu = \sqrt{\mu'\mu}$ and $\sigma$. When given a set of $N$ 2D-points $x = \mu + e$ from a Rice distribution, getting an unbiased estimate of $\sigma$ is straightforward and works as it does for the Rayleigh distribution. The obvious idea to estimate $\nu$ is $\sqrt{\bar{x}'\bar{x}}$, but if I'm not mistaken, the square has expectation $E(\bar{x}'\bar{x}) = \mu'\mu + E(e'e) = \mu'\mu + tr(\frac{1}{N}\Sigma) = \mu'\mu + \frac{2}{N} \sigma^{2}$. And this is before taking the square-root which introduces another bias. The Wikipedia page mentions methods for parameter estimation, but they seem to be for the situation when we don't have the original 2D points, but just their Euclidean norm. It's possible to apply these methods after calculating the distance of the 2D points, they are complicated and based on fewer information than what I have.


title: how to measure the transform ability of a query on Internet Advertising

text: I have a task which aims to measure the transform ablity of a keyword , the transform ablity can be understood as the probability of a query on the search engine brings an order to the customers.For example, I searched "apple" on the search engine then I go the website and buy a notebook, then "apple" bring a transaction to the apple company and the query "apple" has a strong ability to bring a transaction. (Maybe this example is not good but it could explain the idea) Now I have the query log and the click log of advertisement on the search engine , how can I use these data to analyze the transform ablity of a query ? I think maybe some clarification algorithm could give answer but I'm not sure the result though. Could you please give me some hints? Any feedback will be appreciated!


title: Construct artificial slightly overlapping data for PCA plot

text: I am trying to construct artificial data which show two distinct groups in a PCA plot. However, the two groups should still slightly overlap. The following approaches came the closest but I am still not happy with it. Even changing the parameters in `rnorm` and `runif` didn't work out. One group has always a much bigger spread than the other.               d1 <- {}     d2 <- {}     for(i in 1:20) {           d1 <- cbind(d1, rnorm(20, mean=100, sd=10))          d2 <- cbind(d2, rnorm(20, mean=110, sd=5))           # d1 <- rbind(d1, runif(100, 1,4))           # d2 <- rbind(d2, runif(100, 1,3))     }          d <- as.data.frame(rbind(d1,d2))     res <- prcomp(d, center = TRUE, scale = TRUE, na.action = na.omit)     d$group <- as.factor(c(rep("g1", 20), rep("g2", 20)))      ggplot(d, aes(V1,V2)) +                        geom_point(aes(color = group), size=4)


title: Joining results from different laboratories

text: I've done a large number of tests (100). I also have access to a thesis where it is reported the same type of tests on the same type of material. The mean and COV values are reported, along with the histograms. I'd like to use this information, to join it with my data and from there do some model selection. Since the number of tests done by previous researchers (80) is similar to what I did I think Bayesian Theory can only be useful if one of the series of results is considered representative of the (unknown) population. Maybe joining all the data makes more sense. Any advise? Thanks!


title: Simple text classifier: classification taking forever?

text: I work for a small tech startup, and I want to classify or users into demographics based on the domain of their email address. When users sign up to our site, they can enter a job category, or pick "other". The goal is to classify as many of the "other" type as possible using a bag-o-words approach. To do this, I have written some code in Python. For each user, I look at the domain name of their email address and scrape the text from their homepage (using Beautiful Soup). I also look for an "about us" page, which I also scrape. What I'm left with is a map of domains to text. Some domains are classified (i.e., users whose email address comes from this domain have self- classified their job types), and some aren't (those users who have self- classified as "other"). The total data set for classified users is about 2000 (neglecting domains like gmail and hotmail and [I can't believe I'm about to type this] aol). I'm using a train/test split of 75/25. Using `scikits-learn`, I'm trying to implement a simple classifier, but there seems to be an issue with either convergence or performance. The data set doesn't seem particularly big, but the two classifiers I've tried (`Perceptron` and `RidgeClassifier`) seem to be having some issues finding a fit. I haven't really tried to change the parameters for the classifiers, and it's not clear to me which nobs I should be turning. I lack intuition into this problem, and it's difficult for me to tell whether the issues I'm having are due to not enough data, or what. I'd like to know   * Am I barking up the wrong tree? Has anyone tried something like this and made it work?   * Do other ML packages for Python do a better job of text classification? (I'm looking at you, nltk.)   * Is my data set large enough? Are there any "rules-of-thumb" for how much data you'd need for something like this (~5-10 categories)?   * What's a reasonable amount of time for the learning to take? Are there any hints that will tell me the difference between "this is really hard" and "this isn't going to work"?   * I've tried to follow the examples here and here. These examples are pretty speedy, so it makes me worry that I don't have enough data to make things work nicely. Is the "20 newsgroups" classification problem typical, or does it show up because it's easily solvable? Any guidance here would be appreciated! * * * As an update: the huge performance hit seems to come from the "vectorizer": that is, the thing that maps a vector of words to the reals. For some reason, Tfidf was taking a long time to do its thing---I switched to a different vectorizer, and now things run quite quickly. In regards to the actual learning, I've found that the Naive Bayes routines work pretty well out of the box (f-score around 70-75%, which is good enough for now). The model that I found works the best, however, is one based on a linear SVM (`scikits.svm.LinearSVC`), which gets me somewhere in the 80-85% range with a bit of tinkering.


title: How to compute the normalized misclassification error and the misclassification error?

text: How is usually computed the normalized misclassification error ? I would like to know the formula for a 2 class problem of size $s_1, s_2$ and for a k class problem of size $s_1,\dots, s_k$ I would also like to know what formula is used by the r randomForest package to compute the misclassification erron on a $k$ class problem.


title: Fast evaluation of univariate KDE?

text: I have a KDE obtained with a large sample ($n=10000$). Given that this is defined as the mean of the kernel evaluated at $n$ points, the evaluation of this kernel slows down my code (I need to evaluate it thousands of times). I was wondering if there is a trick to avoid taking the mean of the $n$ points every time I evaluate the KDE?


title: What are the degrees of freedom of a distribution?

text: I am dealing right now with a lot of distributions, e.g., $F$, $t$, $\chi^2$. I was wondering why do these degrees of freedom signify for distributions such as the $F(m,n)$ distribution?


title: Zero inflate models vs generalized mixture model

text: Hi I am looking to compare the fit of a zero- inflated mixture model and a poisson mixture model, the random effects in both models are different. Comparing the fitted values of both models ignores the complexity of the models, and using model selection methods such as DIC, AIC etc... is not a straight forward exercise due to the difference in mixture distributions.... so I am wondering if there is a better way? Thank you,


title: Is there any meta-approach for variable selection based of measures of similarity between each two variables?

text: Is there any meta-approach ( or mayby I should say universal approach which works with different measures ) for variable selection which is based on similarity matrix which every entry `S[i,j]` contain value `M(X_i,X_j)` \- mutual information, correlation, reg.coef., R2 of regression etc. between two variables ? In the case of assymetric measures such as KL divergence and conditional enthropy it could be some approach based on algorithms like PageRank.


title: How to cluster temporal trajectories of various social network statistics

text: I have social network data for 100 open source repositories. The data consists of 12 snapshots (1 for every month of a year) of the networks for several repositories. Out of this I want to extract time series for various SNA characteristics/statistics, as in the tables below:               **Table: Repository #1**         Centralization Connectedness Graph Hiearchy     Jan   .5            .8             .1     Feb   .6            .9             .2     Mar   .7            .7             .2          **Table: Repository #2**         Centralization Connectedness Graph Hiearchy     Jan   .3            .5             .1     Feb   .6            .3             .4     Mar   .2            .7             .2          etc...until repository #100      Now I want to identify clusters of projects that have similar social network characteristics, as they progress over time. What is the most suitable technique for clustering projects according to multiple time series simultaneously?


title: Test/score or measure for predicted versus actual long-run proportions/averages?

text: What I am doing is that I have categories A, B, C, D, E. I am computing predicted long-run proportions of falling into each category and computing actual long run proportions of falling into each category. For a vector $longrun = [A, B, C, D, E]$ predicted = $[0.22, 0.33, 0.11, 0.09, 0.25]$ actual = $[0.20, 0.30, 0.10, 0.15, 0.25]$ How do I properly/correctly compare how close these measures came to one another? What score or test/measure should I consider? Many thanks in advance!


title: Testing a dataset against a distribution with parameters estimated from that dataset

text: I am trying to figure out the best distribution to fit some data to, and I'm not sure if what I am doing is statistically correct. My data consists of 20 samples / year over 10 years. For each sample I have run a distribution fitting algorithm (using fitdistr() in R), to get the estimated parameters for each type of distribution. I am testing gamma, chi-squared, weibull and lognormal distributions. My next step was to then run a Kolmogorov Smirnov test, using the sample data, and setting the parameters as estimated from that data. I was going to find which distribution was the overall 'best' (lowest average p-value for all 200 samples), and say that this was the distribution my data described. I have read that using the KS test in this way is incorrect and the resulting p-values will be unreliable. I'm not sure if I can use the KS test in this way, or if I should do and maximum likelihood estimation.


title: Do low $R^2$ values mean that my vector autoregressive model is bad?

text: I have a VAR model, which shows very low $R^2$ values (below 0.05). Does this mean that my model is very bad in explaining? ![enter image description here](http://i.stack.imgur.com/cEDIA.png)


title: Machine learning visualize dataset obtained from machine learning

text: I have analyzed text by finding sentiments of it. Basically have a huge file with text and its associated sentiment. Obtained by doing Machine learning. Used logistic regression algorithm. Now I want to do analysis such as statistics, visualizations, graphs etc. What would be the different kind of analysis I can perform on the data for example, suppose I have twitter data & its sentiments of every tweet then what are the different possibilities of data analysis I can perform? My goal is to convert the dataset to a web app so that users can try it out with interactive visualizations, statistics, graphs etc. Also, it would be helpful if you could suggest even API's or tools.


title: Reported error rates on neural networks

text: It is common to depict the error rates of types of neural networks in a table, for example, see the MNIST website. However, because of the non-determinism caused by weight initialization the actual error rates may vary even under a single setting of hyperparameters. My question is therefore: what statistic over multiple runs is usually depicted? Are the reported error rates the mean error rates? Are they the minimum error rates?


title: Neural Networks and Picture Recognition

text: I have spent a bunch of time looking at this series of videos (Neural Network Tutorial), by Ryan Harris: https://www.youtube.com/watch?v=Q_5B3GuWPCc&index=41&list=PL29C61214F2146796 I am using his libraries to implement my solution. My objective is to be able to train a neural network to detect a picture of a Coca-Cola can. I have several 80 x 108 pixel pictures of the cans. So this means I will need a layer of 6480 inputs. My goal is to have it output to a single output, ideally this output would show on a scale of 0 to 1. The values would show how certain the network is that it is a Coca-Cola can. For example 0.30 would be 30 percent sure. I'm looking for resources or pointers that would show me how this is done. Thanks for any help you can provide


title: Evaluating sum of slopes of a signal

text: I am from India, a novice in Research. Can you just tell me whether I can take a graph of an ECG signal as input and find the slope of this continuous signal at every given point in the graph. I need to get the sum of the slopes from the input ECG signal. Is there tool/method to do this activity of finding slope? I am working in R.


title: rms package, getting zero Dxy in validation of cph model

text: I'm doing Cox proportional hazards analysis using Frank Harrell's rms package (4.2.0), with 3 strata. My model seems reasonably predictive on training data (Dxy=-0.537), but somehow in bootstrap validation Dxy is exactly zero, which doesn't seem likely. I'm specifying `u=10` in `validate.cph` due to the use of strata, but other values of u don't make a difference. Am I doing something wrong?               # some data stats        > summary(d$time)        Min. 1st Qu.  Median    Mean 3rd Qu.    Max.       0.019   9.839  14.780  12.430  14.880  19.940      > with(d, tapply(time, project, summary))     $`1`        Min. 1st Qu.  Median    Mean 3rd Qu.    Max.       0.019   9.725   9.818   8.561   9.875   9.941          $`2`        Min. 1st Qu.  Median    Mean 3rd Qu.    Max.       0.312  10.370  19.770  14.990  19.860  19.940          $`3`        Min. 1st Qu.  Median    Mean 3rd Qu.    Max.        0.191  14.740  14.820  13.380  14.880  14.940     > s <- cph(Surv(time, event) ~ predictor + strat(project), data=d,          x=TRUE, y=TRUE, surv=TRUE)     > s     Cox Proportional Hazards Model          cph(formula = Surv(time, event) ~ predictor + strat(project),         data = d, x = TRUE, y = TRUE, surv = TRUE)                 Status     Stratum     No Event Event       project=1      702   123       project=2      313    78       project=3     2126   190                           Model Tests       Discrimination                                           Indexes     Obs       3532    LR chi2    382.77    R2       0.131      Events     391    d.f.            1    Dxy     -0.537     Center 39.0395    Pr(> chi2) 0.0000    g        1.695                   Score chi2 279.06    gr       5.445                   Pr(> chi2) 0.0000                    Coef   S.E.   Wald Z Pr(>|Z|)     predictor 0.2341 0.0153 15.29  <0.0001          > validate(s, B=20, u=10)                index.orig training   test optimism index.corrected  n     Dxy       0.0000   0.0000 0.0000   0.0000          0.0000 20     R2        0.1306   0.1300 0.1306  -0.0006          0.1312 20     Slope     1.0000   1.0000 1.0056  -0.0056          1.0056 20     D         0.0700   0.0698 0.0700  -0.0002          0.0702 20     U        -0.0004  -0.0004 0.0002  -0.0005          0.0002 20     Q         0.0704   0.0702 0.0698   0.0004          0.0700 20     g         1.6948   1.6939 1.6948  -0.0009          1.6957 20


title: How to make a matrix positive definite?

text: I'm trying to implement an EM algorithm for the following factor analysis model; $$W_j = \mu+B a_j+e_j \quad\text{for}\quad j=1,\ldots,n$$ where $W_j$ is p-dimensional random vector, $a_j$ is a q-dimensional vector of latent variables and $B$ is a pxq matrix of parameters. As a result of other assumptions used for the model, I know that $W_j\sim N(\mu, BB'+D)$ where $D$ is the variance covariance matrix of error terms $e_j$, $D$ = diag($\sigma_1^2$,$\sigma_2^2$,...,$\sigma_p^2$). For the EM algorithm to work, I'm doing dome iterations involving estimation of $B$ and $D$ matrices and during these iterations I'm computing the inverse of $BB'+D$ at each iteration using new estimates of $B$ and $D$. Unfortunately during the course of iterations, $BB'+D$ loses its positive definiteness (but it shouldn't because it is a variance-covariance matrix) and this situation ruins the convergence of the algorithm. My questions are:   1. Does this situation show that there is something wrong with my algorithm since the likelihood should increase at every step of EM?   2. What are the practical ways to make a matrix positive definite? Edit: I'm computing the inverse by using a matrix inversion lemma which states that: $$(BB'+D)^{-1}=D^{-1}-D^{-1}B (I_q+B'D^{-1}B)^{-1} B'D^{-1}$$ where the right side involves only the inverses of $q\times q$ matrices.


title: Relation between electability and probability

text: Is the meaning of 'electability' the same as 'probability' to win? The way I understand it right now, when a survey says "Obama's electability is 63% means: there's X chance (?) that Obama will obtain 63% of the total votes.". [what is that X, by the way, if any]. In that sense electability is not the same as probability (as in the coin flip). Because in the case of coin flip, 63% means: if we perform many many experiments of coin flips, in 63% of occasions Obama will come out as the winner (the number of votes is irrelevant here). Did I get it right? I was trying to get some answers from googling, but so far I haven't found any article that specifically explains the relationship between those two concepts. === UPDATE 1: Mmm.... I think it's the same. Preference by a person can be seen as one experiment (as in coin flip). So..., if they say "63% electability", the meaning is: out of N experiments (persons being surveyed, or coin flip if you will), in 63% of occasions (or 63% of them) Obama is elected (or in the coin- flip equivalent: it is head facing up). Correct? === UPDATE 2: So how exactly do they perform electability survey? I don't think it's as simple as picking, let's say... 2000 persons to be asked, and simply calculates the proportion for each candidate from the answers. In the case of coin flip (bernoulli trial), the bias is already known. For that reason, we can build the (theoritical) probability distribution..., that curve. In the case of electability survey, on the other hand, the bias is unknown. It's the thing we want to find out.... How? I think it's by deduction..., it's deducted from the curve. So, we need to apply some techniques to build that curve, out of sample population. The bias can then be calculated from the properties of that curve. The question now is: what are those techniques to build that curve? ===== Is my line of thinking regarding this subject correct? Thanks in advance, Raka


title: Suppression effect in regression: definition and visual explanation/depiction

text: What is a suppressor variable in multiple regression and what might be the ways to display suppression effect visually (its mechanics or its evidence in results)? I'd like to invite everybody who has a thought, to share.


title: Is there any stochastic method to  approximate the likelihood?

text: I am looking for a list of stochastic algorithms to approximate likelihoods.


title: Significant multiple linear regression model with non-significant betas?

text: In multiple linear regression, let us say that the F test shows that the model is significant. But the t-tests for beta values does not say that the beta values are non-zero. What can we conclude in such a situation? Does the fact that the tests for the beta values failed affect the fact that the model is significant?


title: Sampling dataset, choosing among N dimensions

text: I'm new in stats, and maybe this is a duplicated question, but I could not find a similar one. I'm trying to reduce a dimension of my dataset. Maybe _reduce_ is not a good word. I need to sample some of my dimensions. **Setup** For example:   * I have $M$ events (let's say $M \tilde{} 60$). They are ALL labeled.    * I have $K$ trials/repetitions (let's say $K \tilde{} 10$ or more). $K$ different sets of $M$ events. Therefore, I have a matrix $K\times M$ Now I need to choose $N$ events among $M$ of them as my representatives. I don't need to reduce the dataset, I just need to pick $N$ columns (events) from the matrix $K\times M$. These $N$ events should reflect the behavior of particular trial. I'm interesting in corresponding labels of $N$ chosen events. I know that events are quite different, and probably if there is some correlation it's safe to assume that is a nonlinear and in the best case maybe linear. Further, I want the best representation. Meaning $N$ is not an input (of course $N < M$). **Discussion** I considered PCA, kernel PCA and sampling.   After some reading, I found the following:   * **PCA** will reduce my dataset, but I just need to pinpoint the events. Even if I brake down the algorithm, just choose the appropriate eigenvalues, and reverse back to the original dataset, PCA still remains as linear projection.    * **kernel PCA** is basically nonlinear, but still it will reduce my dataset instead of choosing particular event.    * **sampling** , I went through different types (random, systematic, stratified, cluster ...) but I could not find the particular type. I need something to capture nonlinear connection between the events. Any pointer or explanation of the problem would be much appreciated. Thanks.


title: How to get p-values for model

text: I am using the example from Shumway's book. How do I get the p-values for all the coefficients? Do I divide the coeffients by the standard errors?                Call:         arima(x = cmort, order = c(2, 0, 0), xreg = cbind(trend, temp, temp2, part))              Coefficients:                  ar1     ar2  intercept    trend     temp   temp2    part               0.3848  0.4326  3075.1482  -1.5165  -0.0190  0.0154  0.1545         s.e.  0.0436  0.0400   834.7157   0.4226   0.0495  0.0020  0.0272              sigma^2 estimated as 26.01:  log likelihood = -1549.04,  aic = 3114.07


title: Elasticity of log-log regression

text: In the log-log regression case, $$\log(Y) = B_0 + B_1 \log(X) + U \>,$$ can you show $B_1$ is the elasticity of $Y$ with respect to $X$, i.e., that $E_{yx} = \frac{dY}{dX}(\frac{Y}{X})$?


title: Multiple comparisons in a non-parametric test

text: I am working with a data set that is non-parametric and has 12 treatments. I performed the Kruskal-Wallis test and got a significant $p$-value, and now I would like to conduct a multiple comparisons procedure to see which of the treatments differ significantly. There is a lot of information regarding this topic, but I haven't found anything that specifically addresses this issue. Any ideas??


title: probability of tally having a count of 0 in x variables after sampling with replacement

text: I am a non-statistician testing a randomising algorthm. If I have a sample of size 100 (numbered balls) with replacement and that 50 such samples are independently drawn. To my understanding, the binomial distribution will describe very well the chance that any one particular numbered ball has not been selected in any of the draws (0.605). if after tallying up all the numbers drawn, what is the probability that, say 40 OR MORE of the numbered balls (ANY 40) were not chosen in any of the 50 samples (i.e. any 40 balls have a tally/count of 0)? If possible, could the answer be devoid of R (or other statistical software) arguments, as I wish to understand the logic, rather than just achieve a result. Thank you. ### Edit Thanks to whuber, this question has been altered from its earlier erroneous use of the word "Permutations". In response to the second comment (which will probably not fit as a comment), the exact use of this is as part of a monte-carlo type simulation of a randomising algorithm. I am writing a randomising algorithm in C# which is quite complex, with variable block sizes (as twins of the same gender need to be randomised to the same arm), 2 different interventions and 1 control arm, and a different ratio of the 2 intervention arms. I am running tests on the algorithm before rolling out to trial sites. If a block size is for example 8, with 4 receiving treatment a, 2 receiving treatment b, and 2 control, there will be $$\frac{8!}{4!2!2!}$$ permutations of allocaton order for this block. I wanted to check if any possible permutations were not generated after thousands of runs, and if not, what the probability of having n unused permutations was.


title: Meaning of this expectation equation?

text: I was actually looking at this problem on slide 12. I will write it here briefly: **Problem:** Unknown number of people arriving in a fixed time period and my goal is to maximize my probability of picking the best candidate **Assumptions:**   * Assume people arrive in interval [0,1] independently   * May assume also uniformly   * Can NOT beat probability $\frac{1}{e}$ to win. **Proof:** Fix a wall at time $T$ and select the best candidate after $T$. Let $t$ be the time at which the best candidate arrives, then: $Prob(Win) \geq E_{t}[1_{\\{t>T\\}}\frac{T}{t}] = \int_{T}^{1}\frac{T}{t}dt = -T~ln~T$ This is optimum at $T=\frac{1}{e}$ with a $Prob(Win) = \frac{1}{e} = 36.7\%$ I have a few questions. Can someone please explain how to interpret the expectation equation? In addition, what is the significance of the statement " _Assume people arrive in interval $[0,1]$ independently._ "? Does it mean that I cannot consider some other time interval, say $[0, \alpha]$ or do I need to just multiply $\alpha$ by $\frac{1}{e}$ to get the actual time of making the decision? I am just trying to figure out what I should do to convert this into an actual real-world implementation. Any suggestions?


title: How to compute correlation between/within groups of variables?

text: I have a matrix of 1000 observations and 50 variables each measured on a 5-point scale. These variables are organized into groups, but there aren't an equal number of variables in each group. I'd like to calculate two types of correlations:   1. Correlation within groups of variables (among characteristics): some measure of whether the variables within the group of variables are measuring the same thing.   2. Correlation between groups of variables: some measure, assuming that each group reflects one overall trait, of how each trait (group) is related to every other trait. These characteristics have been previously classified into groups. I'm interested in finding the correlation between the groups - i.e. assuming that the characteristics within in group are measuring the same underlying trait (having completed #1 above - Cronbach's alpha), are the traits themselves related? Does anybody have suggestions for where to start?


title: Creating a cluster map

text: Is it possible to create a cluster map by just using pairwise correlations? If we have $4$ variables than there are $\binom{4}{2}+4$ distinct pairwise correlations. So on both the x-axis and y-axis we would have the four variables. Highly correlated variables would be plotted close together and lower correlated variables would by plotted further apart.


title: Appropriate distributions for infrequent events with copycats

text: I know that infrequent events, especially disasters like earthquakes, can often be modeled by the Poisson distribution. I was thinking of recent events like the Boston bombing and the ricin envelope; these events are infrequent disasters, but aren't independent of past events, because of copycats. What are some good distributions that model infrequent copycat behavior? I imagine that events would come in "packets" with a Poisson distribution, but the size of each packet of events would vary some other way, e.g. a normal distribution. Is there any research on "copycat" disasters?


title: Finding value of residual

text: Consider the LSR line: y = 1158.86 - 5.54x A researcher sought to examine the effect of average teacher salary on the average (total) SAT score. The data, measured at the state level so 50 observations in total, were average teacher salary (measured in thousands of dollars) and average SAT scores out of 1600 points (the old SAT format). Alabama has an average teacher salary (in thousands of dollars) of 31.14 and an average SAT score of 1029. Based on the least squares regression line, what is the value of the residual for the state of Alabama? Can someone explain what the residual is? How do I calculate what it is in this case? Why do they give me information about Alabama specifically (i.e. the salary _and_ score)?


title: vector fit interpretation NMDS

text: So a colleague and myself are using principal component analysis (PCA) or non metric multidimensional scaling (NMDS) to examine how environmental variables influence patterns in benthic community composition. A common method is to fit environmental vectors on to an ordination. The length and direction of the vectors seems somewhat straighforward but I don't understand how an R squared value or a p-value is calculated for these vectors. I have looked at a dozen papers and the most I can gather is that these numbers are calculated using permutations of the data. This does not seem very intuitive. What data is being permuted? How does this permutation create an R squared value and what variance is being explained? My limited understanding of an R squared value comes from linear regressions. I need to explain this to people who have little to no background in statistics so any help understanding these concepts or a link to an available text would be greatly appreciated. Thanks so much!


title: Zero conditional mean assumption

text: How can I check the zero conditional mean assumption for a multiple linear regression? I read that when the assumption is violated, your model is misspecified. And your estimators are biased. So we need to check for it, but I'm completely oblivious as to how this is done.


title: Research paper - section dependent variable for logistic regression

text: In my paper, in the section methodology, I have a subsection called **dependent variable**. But since I am performing a logistic regression, there seem to be 2 (or even 3) things that could be named the dependent variable: First you have your `Y`, which is a measured variable. Second there is the predicted`log-odds of the probability`. Third there is the predicted `probability` (easily calculated from the log- odds). I am wondering which one to describe in my paper in this section? I prefer not to talk about the log-odds, since it makes no sense to make things more complicated at that stage of the paper. I wanted to say the actual Y, cause that's what people are used to in linear regression, but it seems the Y is not _dependent_ of the independent variables, its dependent of the probability. Could I safely talk about the probability, without breaking any rules? Does anyone have any experience with this? Maybe an example of a research paper/article.


title: Backpropagation vs Genetic Algorithm for Neural Network training

text: I've read a few papers discussing pros and cons of each method, some arguing that GA doesn't give any improvement in finding the optimal solution while others show that it is more effective. It seems GA is generally preferred in literature (although mostly people modify it in some way to achieve results they need), then why majority of software solutions seem to use backpropagation only? Is there some general rule of thumb when to use one or another? Maybe it depends on type of NN or there exists some state of the art solution which generally outperforms others? If possible I'm looking for general answers: i.e., "if the NN is huge, GA is better", or "GA is always better but has computational performance issues" etc...


title: Variance for hit-and-miss Monte Carlo method and importance sampling

text: Variance for Hit-and-Miss Monte Carlo is given by $Var(\theta)=\frac{\Theta*(1-\Theta)}{N}$ where $\theta$ is the estimated probability of Hit and N is the number of simulations. Can someone explain why? And what will be the variance when Importance Sampling is used?


title: Identifying instrumental variables for structural model

text: This question is about finding valid instrumental variables. Variables are:   * the response variable for the structural equation, $y$   * the jointly endogenous RHS variable in the structural equation, $m$   * two control variables ($x_1$ and $x_2$) that will always appear on the RHS of the structural equation   * four potential instrumental variables, denoted by ($z_1$, $z_2$, $z_3$, $z_4$) The problem of course, is that you have no idea if the potentional IVs are admissible. Your mission in life is to identifity those Instrumental Variables for which the exclusion restrictions are valid, and to estimate the casual effect of $m$ on $y$ in a linear regression on the form: $$y=x_1\alpha_1+x_2\alpha_2+m\beta+u$$ I dont have any background in these area; is that question trivial or not ?


title: Methodology for validation of stochastic simulations with Kolmogorov-Smirnov test

text: I'm a phd student in Geography, i need some help (or good ressources) to understand why and when i need to use PIT (Probability integral transform) in my validation program for simulation. I explain the context : Actually, i'm running a lot of stochastic simulations (more than 100000) because i use genetic algorithm to optimize parameters values in my model. In this model, we work on system of settlements in old times, and we try to obtain the best hierarchy of settlements. To obtain this result we construct a specific function fitness with multi goal objective. One of this objective is really strict, i take all of my population at the end of the simulation and i try to test the possible **lognormality** of this rank size distribution => result is ok, or not ok... ![Rank size distribution](http://dudarev.github.com/datavis/006_rank_size.png) Actually, i use a classic kolmogorov test to validate/invalidate the lognormality of my city population data, but _ouch_ , this distribution is "unknown" at the end of simulation, so i use this step to resolve this problem : 1 - I try to estimate the standard deviation and the mean of my "unknow distribution" with MLE method (http://en.wikipedia.org/wiki/Maximum_likelihood) 2 - With std and mean for paramter, i generate a new lognormal distribution, named "generated theorical distribution" 3 - I make a uniform transformation on my data, using the "Probability integral transform" or PIT method. _But here i don't understand this step, why i cannot use directly my data, perhaps you can enligth me on this point ?_ ... 4 - I compare the "generated theorical distribution" with the sorted simulation data after the uniform transformation (step3) 5 - After that i extract p-value, and D value to validate / invalidate the lognormality of my simulation data. Now, my question on this method :   * I'm not statistician,so perhaphs these step are false ... in this case, can you explain me the step with the problem ?    * Also, i don't understand the step 3, and the necessity of the "Probability integral transform" of my simulation data here, can you enlight me on this point ?    * Do you think i can use this methodology, and ks.test with only 100 city in my model, so only 100 values in the sample used in lognormality test ? **Update 1** The code which run the PITS :                 def evalKs(sample: Array[Double],dist:LognormalDist)=       {         data = new DoubleArrayList(sample)         v = data.elements         n = data.size         val dataUnif = unifTransform (data, dist)         dataUnif.quickSortFromTo (0, dataUnif.size - 1)         ret = kolmogorovSmirnov (dataUnif)         }      The code for unifTransform in SSJ library GOFStat ( http://www.iro.umontreal.ca/~simardr/ssj/doc/html/umontreal/iro/lecuyer/gof/GofStat.html ) :                  public static DoubleArrayList unifTransform (DoubleArrayList data,                                             DiscreteDistribution dist) {            double[] v = data.elements();            int n = data.size();                 double[] u = new double[n];            for (int i = 0; i < n; i++)               u[i] = dist.cdf ((int)v[i]);            return new DoubleArrayList (u);        }      In the class LognormalDist (by default mu = 0, sigma = 1) :                  public double cdf (double x) {           return cdf (mu, sigma, x);        }             public static double cdf (double mu, double sigma, double x) {           if (sigma <= 0.0)             throw new IllegalArgumentException ("sigma  <= 0");           if (x <= 0.0)              return 0.0;           return NormalDist.cdf01 ((Math.log (x) - mu)/sigma);        }      In the class NormalDist                  public static double cdf01 (double x) {        /*         * Returns P[X < x] for the normal distribution.         * As in J. L. Schonfelder, Math. of Computation, Vol. 32,         * pp 1232--1240, (1978).         */                double t, r;                if (x <= -XBIG)              return 0.0;           if (x >= XBIG)              return 1.0;                x = -x/Num.RAC2;           if (x < 0) {              x = -x;              t = (x - 3.75) / (x + 3.75);              r = 1.0 - 0.5 * Math.exp ( -x * x) * Num.evalCheby (NORMAL2_A, COEFFMAX, t);           } else {              t = (x - 3.75) / (x + 3.75);              r = 0.5 * Math.exp ( -x * x) * Num.evalCheby (NORMAL2_A, COEFFMAX, t);           }           return r;        }      Thanks a lot for your advice, critics on this students work in progress :) SR


title: SVM parameter selection and cross validation

text: Have a quick question about parameter selection for an SVM. I'm using a rbf kernel, so trying to optimize C and gamma. I have an example set of around 4500, about 700 features, and using 700 examples from the set for testing. My dataset does consist of time series. I've been using a 5 fold cross validation with a grid search to find the optimal parameters for the test set and have continued to noticed fairly large differences between the accuracy of my training set vs. the accuracy of my test set. Note, however, when I say accuracy, I have imposed a cost matrix when evaluating the fit of the model such that certain classes have much higher costs when misclassified (note i also ran the svm with unequal class weights). Because my data is a time series, I'm wondering if I should use a different approach from cross validation e.g. a moving window evaluation or something similar. Is cross validation the best approach? Are there other ways to search for the optimal parameters? And also, are there ways to speed up the parameter search (I've heard of using minimum finding algorithms as an alternative to a grid search, which I'm considering implementing)? Any thoughts would be most welcome. Thanks.


title: Marginal Likelihood & Truncated Posterior Mixing

text: It has now become common to use Bayesian inference to find the best solution for exoplanet orbits fitting. In order to find the best solution one has to explore a very large parameter space and some authors propose a technique called _Truncated Posterior Mixture_ to approximate the marginal likelihood. Unfortunately, there seem to be some people acquainted with Bayesian statistics casting some doubts on its efficiency. Could someone acquainted with the subject comment on if the criticism is justified? Many thanks.


title: Statistical test for cumulative curves

text: I have three treatment groups represented in three cumulative curves. What is the best test to show differences among them.Knowing that treatments were observed cumulatively for only seven days.


title: What book about use of machine learning in pharmacy would you recommend?

text: I know that machine learning is a very popular tool in pharmacy. Are there any books that describe use of machine learning in pharmacy?


title: Please guide me to learn SAS

text: I am new to SAS Reporting and want to learn it to use in my job. But could not figure out best book for the same. Please guide me, any book names or tutorials available online so tat I can refer and Learn SAS reporting.


title: Hierarchical ordinal regression (or ranking) with prediction constraints on clusters?

text: I am interested in predicting an ordered outcome of 0,1,2 or 3 (0<1<2<3) for individual responses in a bunch of different clusters. In each cluster $i$ of size $n_i$ there is a single 3, 2 and 1 allocated to group of responses, the remaining $n_i-3$ responses are set to 0. To clarify, the situation is as follows: 3,2,1 votes are given to the 1st, 2nd and 3rd best players in a match and 0 to the remainder, each match represents a cluster. My overall goal is to predict the number of votes obtained by each individual after say, $N$ matches (on a new dataset). Now the first thing I thought to do was a mixed-effects ordinal regression, where each match/cluster has a random intercept. But using this formulation, there is no guarantee that it allocates the 3,2,1 and 0's predictions to individuals in a new cluster in the way that it should. The are two primary ways I thought to solve this problem:   1. Just give the 'expected' number of votes to each player in the new cluster based on the predicted probability of being in each class. (I think this is the way to go)   2. Just allocate the 3,2,1, scores to the individuals in each cluster who have the highest probability of obtaining a 3. The rest get 0's. At this point, my big concern is that what it takes to get votes in one match/cluster could be considerably different from what it takes to get votes in another match, as the 3,2,1 votes must be allocated regardless of the performances of the individuals in the match relative to other matches. I imagine a random intercept for match may alleviate this issue, but is there a better way to deal with it? This is less of a problem with (2) that allocates the 3,2,1 votes regardless, while (1) may allocate less or more than the 6 'expected' votes per match. (But this might not actually be an issue?) Is ordinal regression an appropriate approach? Are the ranking based methods discussed here appropriate? The way I see it is, the method I use doesn't need to allocate the 3,2,1,0 votes, if the method's output is a ranking for each match/cluster. Once the ranking is obtained, it is easy to allocate the votes. Is there anything else I should or need to consider? Cheers.


title: Lagrange Multipliers in practice

text: Say we want to minimize the function $f^2({\bf{x}})$, under the constraint $g({\bf{x}})=0$. * * * The classic solution ( **Method I** ) is to introduce a Lagrange Multiplier, and solve: $$\frac{\partial f^2({\bf{x}})}{\partial x_i}+\lambda \frac{\partial g({\bf{x}})}{\partial x_i}=0$$ One could however use a different method ( **Method II** ), where we redefine the cost function as: $$J\equiv f^2({\bf{x}})+\Lambda g^2({\bf{x}})$$ and then we must solve: $$\frac{\partial f^2({\bf{x}})}{\partial x_i}+\Lambda \frac{\partial g^2({\bf{x}})}{\partial x_i}=0$$ Here $\Lambda$ is a set constant, and is not to be solved for like $\lambda$. Solving leads to a solution that is dependent on $\Lambda$, and is equivalent to method I, when $\Lambda\to\infty$. I think the second method has an advantage when the data is noisy and therefore will not exactly satisfy the constraint, but the choice for $\Lambda$ needs to be made manually. **Does Method II have a name? is it commonly used? what are the main pros and cons, other than the ones mentioned?**


title: MCMC packages in R

text: Is there an R package for MCMC that can   * accept my self-defined (log)likelihood function (can be done in `MCMCpack`)  and   * lets the user define contraints to the proposals (like only accept values from the [0,1] interval, or take only positive values (possible in `rSTAN`) I already definded my model in R code but obviously no MCMC package in R can handle constraints on the parameters like `rstan` and I don't want to recode it in `rstan`.


title: Explanation of the Regression Plot in the Matlab Neural Network Toolbox

text: What does the Regression Plot in the Matlab Neural Network Toolbox show? I thought I understood it when I looked at a univariate regression plot, but I've just plotted one for multivariate regression, and it makes no sense to me. My Neural Network takes in 24 inputs, and gives 3 outputs. The 24 inputs are 24 different angles, and the 3 outputs are the (x,y,z) positions in 3d space. I expect that the function mapping the two is highly nonlinear, but I'm not sure. Here is the regression plot, and here is the performance plot if that helps. What does the regression plot show? I do not understand this graph at all. Surely it cant plot the function because that would require a plot in 27 dimensions (24 inputs + 3 outputs). What do the 4 graphs represent? To me it looks like it is saying that the function is linear, could this be true?


title: Robust Support Vector Regression - robust to outliers

text: I've been reading/looking around for literature on support vector regressions that are relatively robust to outliers. I understand that standard SVRs can be significantly influenced by a few large outliers. From what I've read (and I'm no academic to be sure), there appears to be a number of different approaches. My questions are two: 1) can anyone explain generally how people approach building a robust SVR? And 2) is there any open source code that has implemented some of these methods yet? It would be interesting and educational to dig through them I think. Thanks!


title: Does causation imply correlation?

text: Correlation does not imply causation, as there could be many explanations for the correlation. But does causation imply correlation? Intuitively, I would think that the presence of causation means there is necessarily some correlation. But my intuition has not always served me well in statistics. Does causation imply correlation?


title: LDA projection for classification

text: I am dealing with **2 class LDA classification** problem. During a test phase (after training), I'm trying to project a feature vector to lower dimensional space. **How do we get the projected test feature vector**? Is it   1. **Y = (X-mean) * W**   2. **Y = X * W** Which one of above is true? ( **X** is feature vector, **W** is weight vector obtained during training, **Y** is the resultant projected vector).


title: MaxNWts optimal parameter in nnet R packge

text: I'm fitting a one layer neural network with R. I'm wondering whether there is an optimal approach/practical suggestion to set the MaxNWts parameter. Thanks


title: Correlation between quantitative and qualitative variables

text: I have a dataset composed by 5000 observations. Each observation contains the income per year of a person (from 50 to 50.000.000) and the fact of having a car (yes/no). I would like to check if a correlation exists between these two features. Which test I should run? thanks in advance


title: Criteria for choosing the most appropriate logistic regression model

text: I've fitted 16 logistic regression models to my data and I'm not sure as to which model to choose as my final model. I looked at a couple of things to help me choose my final model. 1) significance of the predictor variables 2) AIC and BIC 3) area under the ROC curve 4) cross validation error. A model with low AIC or BIC value may have a smaller area under the ROC curve than another model with a higher AIC or BIC. And the model with largest area under the ROC curve may have a higher cross-validation error than a model with a smaller area under the ROC curve. I was able to narrow the 16 possible models down to 4-5 models but I can't make the final decision based on the criteria I used (1-4) since they seem to conflict with one another. How can I go about choosing the best logistic regression model for my data?


title: A new piece of clue for document classification?

text: I am working on a document classification problem. I am using the typical vector space model to represent a document as doc-term vector. If document has some term, the vector entry for that term is 1. But with the training data, I can get a piece of information which I think may be useful. It is the term popularity within each category respectively. It is defined as below: $$ \text{Popularity} (\text{term}\ T, \text{category}\ X) = \frac{\text{documents belonging to}\ X\ \text{containing}\ T}{\text{total documents belonging to}\ X} $$ This is calculated with the training set. I believe this value can be interpreted as the **representability** of a term for a certain category, which in other words, is the term's popularity within a certain category. I have a feeling that it could somehow influence the vector weight for a term. But I haven't figured out a reasonable way to use it. If I want to use this popularity value for a new document to be classified, I need to know its category first (or do I?). This is kind of a dilemma. Could anyone shed some light? Could this clue be helpful in the classification? ## ADD 1 Here is my first attempt to use this Popularity value. For each doc to be classified, I make an assumption about its category, and multiply the BIT vector with the term-popularity values for that assumed category. And then, I use some machine learning algorithm to calculate the probability for this doc belong to that assumed category with the new vector. I did this assume-predict for each category. And finally I get a list of probabilities for each category. Then I choose the category with maximum probability as the final prediction. BUT, the result turns out not as good as I expected. It seems even worse than BIT (I am not 100% sure about this.)


title: Obtaining predictions from an mgcv::gam fit that contains a matrix "by" variable to a smooth

text: I just discovered that `mgcv::s()` permits one to supply a matrix to its `by` argument, permitting one to smooth a continuous variable with separate smooths for each of a combination of variables (and their interactions if so desired). However, I'm having trouble getting sensible predictions from such models, for example:               library(mgcv) #for gam     library(ggplot2) #for plotting          #Generate some fake data     set.seed(1) #for replicability of this example     myData = expand.grid(         var1 = c(-1,1)         , var2 = c(-1,1)         , z = -10:10     )     myData$y = rnorm(nrow(myData)) + (myData$z^2 + myData$z*4) * myData$var1 + (3*myData$z^2 + myData$z) * myData$var2 #note additive effects of var1 and var2          #plot the data     ggplot(         data = myData         , mapping = aes(             x = z             , y = y             , colour = factor(var1)             , linetype = factor(var2)         )     )+     geom_line(         alpha = .5     )          #reformat to matrices     zMat = matrix(rep(myData$z,times=2),ncol=2)         xMat = matrix(c(myData$var1,myData$var2),ncol=2)          #get the fit     fit = gam(         formula = myData$y ~ s(zMat,by=xMat,k=5)     )          #get the predictions and plot them     predicted = myData     predicted$value = predict(fit)     ggplot(         data = predicted         , mapping = aes(             x = z             , y = value             , colour = factor(var1)             , linetype = factor(var2)         )     )+     geom_line(         alpha = .5     )      Yields this plot of the input data: ![input data plot](http://i.stack.imgur.com/bE8uC.png) And this obviously awry plot of the predicted values: ![output data plot](http://i.stack.imgur.com/CYr3k.png) Whereas replacing the gam fit above with:               fit = gam(         formula = y ~ s(z,by=var1,k=5) + s(z,by=var2,k=5)         , data = myData     )      but otherwise running the same code yields this reasonable plot of predicted values: ![other output plot](http://i.stack.imgur.com/FvsGv.png) What am I doing wrong here?


title: Cox regression and time scale

text: Does X (hazard) variable in Cox proportional hazard regression analysis always have to be time? If not, could you provide an example, please? Can age of cancer patient be a hazard variable? If so, can it be interpreted as the risk of getting cancer at a certain age? Would Cox regression be a legitimate analysis to study the association between gene expression and age?


title: Probability with fair dice

text: If you roll 5 standard 6 sided dies, what's the probability that you will get at least three 2s? I imagine it would be 1 - P(0 twos) - P(1 two) - P(2 twos), but I don't know how to calculate the probability of these.


title: What is the best statistical test for a time series?

text: I have a simple time series with 5-10 data points per data set at regular intervals. I am wondering what is the best way to determine whether two data sets are different. Should i try t-tests on each data point, or look at the area under the curves or is there some kind of multivariate model that would work better?


title: Do you know of any public circular/angular dataset?

text: I've been trying to locate and access circular/angular datasets for my research. I'm particularly interested in datasets that distribute according to a mean centered circular distribution with proper references for citation. Also, some insights in finding this type of data for the future will be appreciated.


title: Naive Bayes model

text: i am trying to implement the naive bayes algorithm in java. But i have a few confusions on the naive bayes model. What is the model of naive bayes? Is it a table ?how we can make prediction from the model? I can implement the algorithm in java code and make prediction.But i heard that's not the right way to do it. we should make prediction from the model ..but how ??If anyone have some idea about this,kindly please try to give an answer to this . And also, how can we split the training set and test set from the dataset ??? Thank you.


title: How to identify patterns/associations in network traffic?

text: The web-based company I work for has a system that receives packets of information from external sites. I would like to know what approaches can be used to identify patterns or associations between the external sites. For example site A sends a packet on subject 1. Shortly after many other sites may send through packets on the same ID. This may be because Site 2 also received the information from Site 1 and then sent it through to us. To make it more complicated, it is also possible that once we receive the data we would send it to one or more sites, which may then re-start the process. The packets contain an ID field that can be used to uniquely identify the subject each packet refers to. There are hundreds of sites, many of which aggregate packets from each other before sending them on to us. Many of the smaller sites could be excluded from the analysis if that simplifies it. We have a lot of data. Over the course of a few days we may see hundreds of packets all relating to the same ID. Also, what approaches could be used to best visualise or display the results/patterns to a non-statistical audience? Order _may_ be important. For example if Site 1 always sends through an ID before Site 2, then Site 2 is (possibly) getting its packets from Site 1. Apologies for the vagueness of the question. I will aim to clarify this question if more details are requested. Thanks EDIT: It would also be of interest to look at the analyses for individual sites.


title: Fisher information metric for hierarchical Bayesian model is negative-definite?

text: I'm strungling at the computation of Fisher information matrix for the hierarchical Bayesian model. For simplicity, consider theta following hierarchical Bayesian model: $$X|\sigma \sim N(0,\sigma^2)$$ $$\sigma|\gamma\sim LN(\gamma,0)$$ $$\pi(\gamma)\propto1$$ $X$ is some observable data point, $\pi(\gamma)\propto1$ means that unknown parameter $\gamma$ has flat prior distribution, and $LN$ denotes the log-normal distribution. The likelihood of the data then is as follows: $$f(X|\sigma,\gamma)=\frac{1}{\sigma}e^{-\frac{x^2}{2\sigma^2}}\frac{1}{\sigma}e^{-\frac{(ln\sigma-\gamma)^2}{2}}$$ Log-likelihood is then: $$L(X|\sigma,\gamma)=-2ln\sigma-\frac{x^2}{2\sigma^2}-\frac{(ln\sigma-\gamma)^2}{2}$$ First order derivatives: $$\partial_{\sigma}L=-\frac{2}{\sigma}+\frac{x^2}{\sigma^3}-\frac{1}{\sigma}(ln\sigma-\gamma) $$ $$\partial_{\gamma}L=ln\sigma-\gamma.$$ Second order derivatives are then: $$\partial_{\sigma}^2L=\frac{1}{\sigma^2}-\frac{3x^2}{\sigma^4}+\frac{1}{\sigma^2}(ln\sigma-\gamma) $$ $$\partial_{\sigma,\gamma}^2L=\frac{1}{\sigma}$$ $$\partial_{\gamma}^2L=-1$$ Thus, the hessian of the log-likelihood is: $$Hess(L)=\begin{pmatrix} \frac{1}{\sigma^2}-\frac{3x^2}{\sigma^4}+\frac{1}{\sigma^2}(ln\sigma-\gamma)&\frac{1}{\sigma} \\\ \frac{1}{\sigma}&-1 \end{pmatrix}$$ The Fisher information metric is defined as follows: $$G=-E_{p(x,\theta)}\left [ Hess(L)) \right ]$$ Since $E_{f(x,\theta)}\left [x^2\right ]=\sigma^2$ we have that Fisher metric is $$G=\begin{pmatrix} -\frac{-2+ln\sigma-\gamma}{\sigma^2}&-\frac{1}{\sigma} \\\ -\frac{1}{\sigma}&1 \end{pmatrix}$$ But (and there comes the twist), if we plug-in parameter values $\sigma=1$ and $\gamma=-3$, we obtain negative-definite matrix: $$G=\begin{pmatrix} -1&-1 \\\ -1&1 \end{pmatrix}$$ But Fisher information matrix has to be semi-definite. So, my question is - what am I doing wrong? I have a hunch that expectation of the parameter $\sigma$ has to be taken having in mind that it is also a random variable, not just the data $X$, i.e. for example $E_{f(X,\sigma,\gamma)}\left [\sigma \right ]\neq\sigma$ but $E_{f(X,\sigma,\gamma)}\left [\sigma \right ]=e^{\gamma+\frac{1}{2}}$ - as $\sigma$ has log-normal distribution. I search for literature on Fisher information matrix formation for hierarchical models but in vain.


title: Adding exogeneous variables to a GARCH model

text: **I am somewhat new to R and am currently stuck with the following problem:** I have two variables, $X$ and $Y$. In the first step I do a GARCH(1,1) fit on $Y$. In the second step I would like to input $X$ as an exogenous variable into the model. Afterwards it would be useful if I could compare both models (the goodness of the respective fit) by AIC and BIC, and if possible save the residuals for both models. As far as I know `rugarch` would be the correct package to use. Could anybody perhaps be so kind to help me with the code - please :) This way I would have a working example and could build on it and learn from it.


title: Is there a smooth probability density function with finite moments and a closed form quantile function?

text: I am looking for a smooth probability density function with finite moments and closed form quantile function. As one knows, an example of smooth probability density function with finite moments is the Gaussian density funciton, but it does not have closed form quantile function. And an example of smooth probability density with closed form quantile function is the Cauchy density function, but it does not has finite moments. Finally, an example of probability density function with finite moments and closed form quantile function is derived by extending the exponential distribution density function by flipping it along the y axis and dividing it by 2, but the derived density function is not smooth at one point (the mean). I wonder whether one can construct a probability density function which is smooth, has finite moments and has closed-form quantile function. Edit: where 'smooth' means _infinitely differentiable_ and 'has finite moments' means _all positive integer moments are finite_


title: Standard Error of Standard Deviations Estimated with gls in R

text: In the gls fit shown below, the estimates of the standard deviation for each level of X are apparently given by the product of (1.000000, 3.913972, 10.684698, 11.350910, 26.476561, 27.255072) times the residual standard error of 0.04896334. How does one estimate the standard error of these estimated standard deviations (or variances)?               > m         X     Y  F     1   1  1.07  1     2   1  1.01  1     3   1  0.99  1     4   1  1.09  1     5   1  0.94  1     6   1  1.00  1     7   1  1.01  1     8   1  0.98  1     9   1  1.00  1     10  1  1.03  1     11  4  3.66  4     12  4  3.75  4     13  4  3.77  4     14  4  3.92  4     15  4  4.08  4     16  4  3.99  4     17  4  3.95  4     18  4  4.10  4     19  4  3.88  4     20  4  4.04  4     21 10 10.13 10     22 10 10.20 10     23 10  9.77 10     24 10 10.28 10     25 10  8.71 10     26 10  9.79 10     27 10  9.82 10     28 10  9.85 10     29 10 10.07 10     30 10  9.63 10     31 20 20.22 20     32 20 19.46 20     33 20 19.02 20     34 20 20.06 20     35 20 20.94 20     36 20 19.92 20     37 20 19.96 20     38 20 20.04 20     39 20 19.67 20     40 20 19.96 20     41 30 31.04 30     42 30 31.40 30     43 30 31.84 30     44 30 30.77 30     45 30 32.13 30     46 30 31.17 30     47 30 30.36 30     48 30 29.95 30     49 30 30.74 30     50 30 30.67 30     51 40 41.14 40     52 40 40.29 40     53 40 42.77 40     54 40 38.36 40     55 40 39.17 40     56 40 39.61 40     57 40 40.73 40     58 40 39.42 40     59 40 40.72 40     60 40 40.24 40     > Fit.gls <- gls(Y ~ X,weights=varIdent(form = ~ 1 | F),data=m)     > summary(Fit.gls)     Generalized least squares fit by REML       Model: Y ~ X        Data: m             AIC      BIC    logLik       78.96207 95.44562 -31.48104          Variance function:      Structure: Different standard deviations per stratum      Formula: ~1 | F       Parameter estimates:             1         4        10        20        30        40       1.000000  3.913972 10.684698 11.350910 26.476561 27.255072


title: getting rid of the batch effects

text: I have a two data set one with $n=15$ another with $n=25$. Each set of data has around 90 variables. The data collected was experimental so when I try to combine my two datasets I see clear differences in some variables due to this batch effect. I see these variables in a scatterplot matrix and some of them do not show separation based upon batches while others clearly show a separation. Is there a way to normalize them to combine these two sets. I tried taking z scores on each of the set separately and combine them based upon their score. When I try to see scatterplot matrix of the z-scores, the batch effects are not visible. Is it a valid way to normalize?


title: How does LASSO select among collinear predictors?

text: I'm looking for an intuitive answer why a GLM LASSO model selects a specific predictor out of a group of highly correlated ones, and why it does so differently then the best subset feature selection. From the geometry of the LASSO shown in Fig 2 in Tibshirani 1996 I'm led to believe that LASSO selects the predictor with the greater variance. Now suppose that I use best subset selection with 10 fold CV, to obtain 2 predictors for a logistic regression model and I have reasonable prior knowledge that these 2 predictors are optimal (in 0-1 loss sense). The LASSO solution favors a less parsimonious (5 predictors) solution with greater prediction error. Intuitively, what causes the difference to arise? Is it because of the way LASSO selects among correlated predictors?


title: Which distribution to use to model web page read time?

text: I have a function that returns an average wait time for a web user. That is, it gives an average time that an average user may stay on a web page, given the web resource length in words. I want to use this function (and the resulting average) in conjunction with a distribution to model an 'average web user' browsing the web. Which distribution(s) might be appropriate for this, and why? Edit: I would also like to particularly know the viability of using an exponential distribution for this purpose. Thank you


title: Market participation of small holder famers

text: I am trying to estimate selection and outcome equations using Heckman two step selection model on Stata. For selection equation the dependent variable is participation (1,0) and independent variables are age, sex, household size, education level, quantity produced. For outcome equation dependent variable is fraction marketed and independent variables are age, sex, household size, education level, quantity produced, distance to nearby market, unit price, marketing experience and household income. How can I exactly type Stata command and run the analysis?


title: Explaining MCMC sampling for a Multinomial Distribution and missing at random

text: I understand how MCMC works, and I understand how Multinomial Distribution works. I have a dataset some of the data are missing at random (MAR). I cannot connect these two dots together (MCMC -> Multinomial Distribution) to draw random samples. Can somebody please explain how MCMC sampling works for a Mutlinomaial distribution using a simple example please. Thanks for your time. Hope I am not asking something very complicated.


title: Reparameterizing the multivariate normal model in Bayesian atatistics, Jacobian of the covariance variance matrix?

text: For the multivariate normal model in Bayesian statistics, where the sampling distribution is $y|\theta,\Sigma \sim N_n(\theta, \Sigma)$, and $\theta$ is a n-dimensional vector, and $\Sigma$ a n by n covariance-variance matrix, I'm trying to re-parametrize the $\Sigma=\Psi^{-1}$. I'm having troubles with finding the Jacobian of $\Psi^{-1} $. How do I find it and what is the answer?


title: How to test the influence of external factors on time series

text: I need an **advice/ best method** on how to rigorously test the **influence of categorical predictors** (states), or their duration (duration of war/peace in this case) **on time series**. **How to prove that the increase in prices of oil is caused by war?** I would like to get some _p-value_.               # Example - influence of war on prices of oil     # (completely speculative of course ;D)     price <- c(13,13.5,13,13.5,13,15,15.5,15.5,15,12.5,12,12,12,12.5,14,15,                14,13.5,12.5,13,13,14,15,13.5,12,12,12,12.5,12,12)     year <- seq(from = 1960, length = 30)     era <- rep(c("peace","war","peace","war","peace","war","peace"),                  times = c(5,4,5,2,3,4,7))     my.data <- data.frame(price = price, year = year, era = era)          # plot     plot(my.data$year, my.data$price, xlab = "year", ylab = "price of oil",          pch = 20, ylim = c(11,17), cex = 2.2)     lines(my.data$year, my.data$price, col = "blue", lty = 2, lwd = 2)      Can anybody help me? Thank you for any guidance in advance! At least the name of appropriate technique will be helpful. * * * ![enter image description here](http://i.stack.imgur.com/Y4A5A.png)


title: Bonferroni-correction after one-way ANOVA

text: I conducted a one-way ANOVA: I had two groups and tested, if they differ in three different DVs. Since I tested three times, I thought it was wise to go for Bonferroni- correction. But then I got the notification that the correction can't be done because there are less than three groups. Can someone explain this to me? I thought the Bonferroni-correction was exactly for this case: testing data on the same groups several times.


title: Modeling financial time series

text: Price time series are not stationary. So we difference them and get the return time series, which are stationary. Does this mean, it is always a good idea to model only the return series of financial assets. Alternatively, do we not need to model prices ever?


title: Scatter plots in ggplot2 with non linear trend line

text: I'm trying to create a scatter plot using ggplot. The x value goes roughly from 0-20, and the y value goes from 0-~50. The data is raw so it's not nice and neat and I need to plot a curved line through the majority of the data. The code I'm trying to use is               ggplot(temperature, Number.of.Fish,aes(x=Temperature,y=Number.of.Fish))+geom_smooth      But I keep getting an error message: error in inherits (mapping, "uneval"), could not find function "Number.of.Fish"


title: Conditional logit model and price elasticities

text: I have a dataset with 20 products and whether someone bought them, their price and other attributes. I am trying to find the own and cross-price elasticities of these 20 goods. I have used a conditional logistic model in Stata and now need to find the elasticities preferably in a matrix. I am thinking of using `margins` or `mfx`. Using `mfx` would I need to use the below? Or is there a better approach to finding the price elasticity matrix in Stata?               mfx predict(bought = 1) dyex(ln(product = '1'))     mfx predict(bought = 2) dyex(ln(product = '2'))     etc.      Regression results were generated as per:                clogit bought productX attrY for X = 1..20 and Y = 1..15.      Could I use `mfx eyex` as to get the own and cross-price elasticities?


title: High p-values in Probit model

text: Can someone help me interpret those results? I get very high p-value in my probit model but I do not understand why... ![enter image description here](http://i.stack.imgur.com/S0wAw.png) Some details about my variables: confucius = dependent (binary 1=yes, 0=no) totalagri = continuous (tons) bribes = percentage size = continuous (km2) coast = binary agrigdp = percentage impchina, expchina, gdp = continuous (value in $) lifeexp = continuous (years) urbpop = percentage coalprod, nargasprod, oilreserv = continuous diamond = binary cropprod = index Thanks!


title: Statistical comparison of data dispersion between two groups

text: I was trying to find a suitable way to compare the relative data dispersion between two groups and also to investigate if they are statistically significant. During the last weeks, I have searched a lot and read several papers and finally come up with a solution. As I am not from statistical background, I would like your advice on this matter, if I am doing all wrong!   1. According to Jeerapa [1], there are several ways to compute confidence interval for a coefficient of variation ($cv=\frac{\sigma}{\mu}$, where, $\sigma$ represents standard deviation, and $\mu$ represents mean value). I have used both Miller's Confidence interval and McKay's confidence interval. The equations can be found in the cited paper (3 and 4).   2. According to Altman [2], we can compute both p-value (significantly different than 0) and standard error from a given 95% CI and mean value.    3. According to this website, we can compute t-statistics of mean differences and then can obtain p-value from the respective t-distribution table. My questions are:   1. Is it a valid way to compare the relative data dispersion between two groups? Can we use Bonferroni correction for more than two groups?   2. Basically, I have normalised a dataset using different methods and compared the relative dispersion. For example, a raw feature's cv is 12%, after normalisation this feature using method a1, the cv became 10%, and after normalisation it using method a2, the cv became 6%. As the CV is not computed from independent dataset, can we perform t test instead of paired t test? For paired t test, we need standard deviation of the differences which we cant compute (cv is a scalar value!). Using the above way, I found the p-value is less than 0.05 for method a1 and a2 where a2's cv is significantly lower than a1's cv.   3. Limitations of CV are: a) same std but different mean would create different cv-value although both data dispersion are same. b) very sensitive when mean close to 0. Most of the average values from method a1 are close to 0 (e.g., 0.6). Is there any alternate way to compare the data dispersion between two or more groups? I would really appreciate your help. [1] Sappakitkamjorn, J. and Niwitpong, S. 2013, 'Confidence intervals for the coefficients of variation with bounded parameters', International Journal of Mathematical, Computational Science and Engineering, Vol. 7(9). [2] Altman, D.G. and Bland, J.M. 2011, 'How to obtain the P value from a confidence interval', BMJ, 2011.


title: Experiment with two groups, pre- and post- treatment assessments

text: I am conducting an experiment with a control group and a treatment group. Both groups will take a pre and post treatment depression inventory. What statistical procedure do I use to determine if the treatment reduced depressive symptoms significantly? T-test or ANOVA, which kind and why?


title: Confidence interval and probability - where is the error in this statement?

text: If someone makes a statement like below: > "Overall, nonsmokers exposed to environmental smoke had a relative risk of > coronary heart disease of 1.25 (95 percent confidence interval, 1.17 to > 1.32) as compared with nonsmokers not exposed to smoke." > > What is the relative risk for the population as a whole? How many things are > connected with coronary heart disease? Of the vast number of things that can > be tested, very few actually are connected to coronary heart disease, so the > chance that any particular thing chosen at random is connected is > vanishingly small. Thus we can say that the relative risk for the population > is 1. But the quoted interval does not contain the value 1. So either there > actually is a connection between the two things, the probability of which is > vanishingly small, or this is one of the 5% of intervals that do not contain > the parameter. As the latter is far more likely than the former it is what > we should assume. Therefore, the appropriate conclusion is that the data set > was almost certainly atypical of the population, and thus no connection can > be implied. > > Of course, if there is some basis for assuming that more than 5% of things > are linked to coronary heart disease then there might be some evidence in > the statistic to support the suggestion that environmental smoke is one of > them. Common sense suggests that this is unlikely. What is the error in their reasoning (as all health organizations agree that there is significant literature regarding the damaging effects of second-hand smoking)? Is it because of their premise that "Of the vast number of things that can be tested, very few actually are connected to coronary heart disease"? This sentence may be true for any randomly chosen factor (ie. how many dogs a person owns with the risk of coronary artery disease) but the a priori probability is much higher for second hand smoking and coronary heart disease than just 'any random factor'. Is this the correct reasoning? Or is there something else?


title: What are key differences between homoscedasticity and homogeneity?

text: Homoscedasticity and homogeneity assumptions are popular and perhaps deal with ANOVA and regression. These assumptions create lot of confusion.


title: Is it 'legal' to publish a meta-analysis online?

text: I would like to present a meta-analysis on a homepage. This should include a forest plot to visualize the effect sizes from individual studies as well as the summary statistics. I was wondering whether this plot is 'legal' as I am taking results (effects sizes) from published papers and present them 'for free' as part of this meta-analysis?


title: Differences between SVR with a linear kernel and linear least squares

text: I've been working on a toy problem of predicting reviews a product will get in the future. I found that SVR with a linear kernel worked better than doing a linear least squares regression on the data (even though both were pretty terrible), but I had no intuition about which one works better and why. What are some examples of using one versus the other? Can both usually be used?


title: Symmetric Kullback-Leibler divergence OR Mutual Information as a metric of distance between two distributions?

text: I need some metric of divergence of two distributions. (They are complex and don't fit with exponential family, normal, log-normal, power-law. Maybe some mixture of that, but I'm not feeling right now figuring that out.) I'm thinking between Kullback-Leibler divergence and Mutual Information. I don't have any default distribution, therefore I don't like DKL to be asymmetric (Of course I can symmetrize it, but I don't feel confident about it.) Which one would you choose? Any other ideas are also welcome!


title: Error in adding interaction in Cox model?

text: I'm doing a survival analysis and after plotting the Schoenfeld residuals and test the significance of the correlation of residuals with time, I've decided to incorporate a interaction in the model. It makes sense since the variable outcome have a varying effect over time, for those who defaulted the previous episode, the probability to have another episode is much quicker than for those who completed the previous episode. I've tried to add an interaction but this is a bit confusing and leads to a **numerical singularity**. > modf <\- coxph(Surv(time, status) ~ Vac + CliForm + Sit + Prison + HIV + age > + Sit:time, data2) > > Warning message: In fitter(X, Y, strats, offset, init, control, weights = > weights, : Ran out of iterations and did not converge However using this also gave the same error: > modf <\- coxph(Surv(ini, end, status) ~ Vac + CliForm + Sit + Prison + HIV + > age + Sit:end, data2) > > Warning message: In fitter(X, Y, strats, offset, init, control, weights = > weights, : Ran out of iterations and did not converge Another option is to separate the time in two:               data4$statusCured <- ifelse(data4$time >= 250, data4$status, 0)         data4$statusDef <- ifelse(data4$time < 250, data4$status, 0)     summary(coxph(Surv(time, statusDef) ~ Vac + CliForm + Sit + Prison + HIV + age, data4))     summary(coxph(Surv(time, statusCured) ~ Vac + CliForm + Sit + Prison + HIV + age, data4))      **My main question is what is the best approach to this problem?** How should I include an interaction in my model? In the case that the variable outcome has a variation with time. Any pointers would be extremely helpful.


title: How to transform observed variables to their underlying latent variable in SPSS?

text: My questionnaire includes 48 questions (observed variables) that represent 8 different factors (latent variables). All the variables are continuous. I need to compute the latent variables before doing the correlation & regression analyses but the problem is that SPSS doesn't directly create latent variables like SEM softwares. How can I create latent variables in SPSS?


title: What is class correlation?

text: I am referring to `iris.arff` dataset that comes with Weka distribution. There are 3 classes and 4 attributes. It gives below information regarding the 4 attributes-               %                  Min  Max   Mean    SD   Class Correlation     %    sepal length: 4.3  7.9   5.84  0.83    0.7826     %     sepal width: 2.0  4.4   3.05  0.43   -0.4194     %    petal length: 1.0  6.9   3.76  1.76    0.9490  (high!)     %     petal width: 0.1  2.5   1.20  0.76    0.9565  (high!)      How is the last column computed? After googling class correlation I found Intraclass_correlation Not sure if it's the same thing, also how can a single class correlation number be given when the number of classes is three?


title: Presentation of summary log-transformed data aiming at easier interpretation

text: I have a list of measurements that is not normal and need to be log- transformed so that parametric statistics can be applied (e.g. mean, CI). I can calculate the mean and CI without problem when the data is transformed, but the way of interpreting the log-transformed value is not straightforward to my readers. (the CI for the object length is 1-5 log cm, people usually will ask what is log cm?) I am tempted to revert the CI with transformed unit back to it original unit, but I wonder if this is the right way to do so, or if there can be other alternatives? The R script below will explain my question fully:               library(epicalc)          # creation of sample data     set.seed(1)     sample.measurement <- 10^(rnorm(1000, mean=1, 0.5))          # plotting histogram for sample data, too skewed to use mean to summarize data,      # need to log-transform it     hist(sample.measurement)     # plotting histogram for sample data, much better this time     hist(log(sample.measurement, base=10))          # the calculation below should be wrong, as it is not normal data     ci.numeric(sample.measurement)          # this one should be better, as the data is log-transformed and become normal,      # but the value is difficult to interpret/comprehend     ci.numeric(log(sample.measurement, base=10))          # I wonder if the method below is appropriate, for the sake of easier interpretation      # using the original units     10^ci.numeric(log(sample.measurement, base=10))