title: Motivation for Kolmogorov distance between distributions

text: There are many ways to measure how similar two probability distributions are. Among methods which are popular (in different circles) are:   1. the Kolmogorov distance: the sup-distance between the distribution functions;   2. the Kantorovich-Rubinstein distance: the maximum difference between the expectations w.r.t. the two distributions of functions with Lipschitz constant $1$, which also turns out to be the $L^1$ distance between the distribution functions;   3. the bounded-Lipschitz distance: like the K-R distance but the functions are also required to have absolute value at most $1$. These have different advantages and disadvantages. Only convergence in the sense of 3. actually corresponds precisely to convergence in distribution; convergence in the sense of 1. or 2. is slightly stronger in general. (In particular, if $X_n=\frac{1}{n}$ with probability $1$, then $X_n$ converges to $0$ in distribution, but not in the Kolmogorov distance. However, if the limit distribution is continuous then this pathology doesn't occur.) From the perspective of elementary probability or measure theory, 1. is very natural because it compares the probabilities of being in some set. A more sophisticated probabilistic perspective, on the other hand, tends to focus more on expectations than probabilities. Also, from the perspective of functional analysis, distances like 2. or 3. based on duality with some function space are very appealing, because there is a large set of mathematical tools for working with such things. However, my impression (correct me if I'm wrong!) is that in statistics, the Kolmogorov distance is the usually preferred way of measuring similarity of distributions. I can guess one reason: if one of the distributions is discrete with finite support -- in particular, if it is the distribution of some real- world data -- then the Kolmogorov distance to a model distribution is easy to compute. (The K-R distance would be slightly harder to compute, and the B-L distance would probably be impossible in practical terms.) So my question (finally) is, are there other reasons, either practical or theoretical, to favor the Kolmogorov distance (or some other distance) for statistical purposes?


title: Computing mutual information

text: I have a problem when computing the mutual information between two variables. Let's consider the following table:                            k=1       k=2       k=3         c = 1         10        20         5         c = 2          5         3        20      I want to calculate the mutual information between the class ($c$) and the cluster ($k$). E.g., for $\text{MI}(c=1, k=1)$, I will calculate the $P(c,k)$ as follows:                            k=1            k~=1         c=1        10/63           25/63         c~=1        5/63           23/63      and I can calculate the MI based on the following: $\text{MI}(c,k) = \sum_{c\in\{0,1\}} \sum_{k\in\{0,1\}} p(c,k)log_2 \frac{p(c,k)}{p(c)p(k)}$ In this case: I am getting the same MI for $(c=1, k=1)$ and $(c=2, k=1)$. What is the reason for this? How can I calculate the MI for each $(c,k)$ based on the first table?


title: Interpretation of a linear regression

text: I have some problems in understanding the results of two linear regressions on the same data. I have 2 variables and I want to know if there is a linear correlation between the two. Let's say velocity of a car and perceived velocity of the conductor. Let's say I have 4 velocities and 2 evaluations of the perceived velocities performed in 2 trials by 15 participants, on a 7 points Likert scale where 1 is very slow and 7 = very fast. Now if I take the AVERAGE values of the evaluations given by the 15 participants for the 4 velocities I get a clear linear relationship with $R^2$ = 0.7. However if I perform the analysis on the raw data (i.e. 15 participants * 2 trials = 30 values for each velocity) I get a $R^2$ = 0.2. Plotting the graph of the velocity data vs perceived velocity I get a mess of points very widely distributed. Why do I get this behavior? Which of the two analyses should I use? The first or the second? How can I interpret my values? Please let me know Thanks in advance


title: Response distribution dependent bias in random forest regression

text: I am using the randomForest package in R (R version 2.13.1, randomForest version 4.6-2) for regression and noticed a significant bias in my results: the prediction error is dependent on the value of the response variable. High values are under predicted and low values are over predicted. At first I suspected this was a consequence of my data but the following simple example suggests that this is inherent to the random forest algorithm:               n = 1000;      x1 = rnorm(n, mean = 0, sd = 1)     response = x1     predictors = data.frame(x1=x1)      rf = randomForest(x=predictors, y=response)     error = response-predict(rf, predictors)     plot(x1, error)      I suspect the bias is dependent on the distribution of the response, for example, if x1 is uniformly distributed, there is no bias; if x1 is exponentially distributed, the bias is one-sided. Essentially, the values of the response at the tails of a normal distribution are outliers. It is no surprise that a model would have difficulty predicting outliers. In the case of randomForest, a response value of extreme magnitude from the tail of a distribution is less likely to end up in a terminal leaf and its effect will be washed out in the ensemble average. Note that I tried to capture this effect in a previous example, "RandomForest in R linear regression tails mtry". This was a bad example. If the bias in the above example is truly inherent to the algorithm, it follows that a bias correction could be formulated given the response distribution one is trying to predict, resulting in more accurate predictions. Are tree based methods, such as random forest, subject to response distribution bias? If so, is this previously known to the statistics community and how is it usually corrected (e.g. a second model that uses the residuals of the biased model as input)? Correction of a response dependent bias is difficult because, by nature, the response is not known. Unfortunately, the estimate/predicted response does not often share the same relationship to the bias.


title: Why isn't RANSAC most widely used in statistics?

text: Coming from the field of computer vision, I've often used the RANSAC (Random Sample Consensus) method for fitting models to data with lots of outliers. However, I've never seen it used by statisticians, and I've always been under the impression that it wasn't considered a "statistically-sound" method. Why is that so? It is random in nature, which makes it harder to analyze, but so are bootstrapping methods. Or is simply a case of academic silos not talking to one another?


title: General question on oversampling

text: I once read (unfortunately I cannot find the source anymore) that oversampling is only "allowed" when the class distribution is very, very skewed, i.e. 1:99 or worse. So in my case I have a class distribution of 30:70. When don't apply oversampling in my neural network (10-20-1, sigmoid (hidden layer), softmax(output layer)) then the f-measure is around 0.10 and the kappa around 0.1, which means in fact nothing to predict/features very bad. When I apply oversampling, so the class distribution is almost 50:50 the f-measure goes up to 0.9 and the kappa to 0.85. As this sounds too good to be true, I think oversampling in this case isn't valid. What is your opinion this? And do you have a good reference for rules of thumb on oversampling?


title: Which model should I prefer for time series forecasting?

text: I have time series as               0.4385487 0.7024281 0.9381081 0.8235792 0.7779642 1.1670665 1.1958634 1.1958634 0.8235792 0.8530141 0.8802216 1.1958634 1.1235897 1.3542734 1.3245534 0.9381081 1.1670665 1.1958634 0.8802216 1.3542734 1.1670665 4.9167998 0.9651803 0.8221709 1.1070461 1.2006974 1.3542734 0.9651803 0.9381081 0.9651803 0.8854192 1.3245534 1.1235897 1.2006974 1.1958634 0.4385487 1.3245534 4.9167998 1.2277843 0.8530141 1.0018480 0.3588158 0.8530141 0.8867365 1.3542734 1.1958634 1.1958634 0.9651803 0.8802216 0.8235792 4.9167998 1.1958634 0.9651803 0.8854192 0.8854192 1.2006974 0.8867365 0.9381081 0.8235792 0.9651803 0.4385487 0.9936722 0.8821301 1.3542734 1.1235897 1.6132899 1.3245534 1.3542734 0.8132233 0.8530141 1.1958634 1.2279813 0.8354292 1.3578511 1.1070461 0.8530141 0.9670581 1.1958634 0.7779642 1.2006974 1.1958634 0.8235792 1.3245534 0.5119648 2.3386331 0.8890464 0.8867365 4.9167998 1.2006974 1.2006974 0.6715839 4.9167998 0.7747481 4.9167998 0.8867365 1.2277843 0.8890464 1.2277843 0.8890464 1.0541099 0.8821301       I am using package "itsmr"-autofit(),"forecast"-auto.arima(),"package"-- functions   1. Autoregressive model              > ar(t)          Call:         ar(x = t)              Order selected 0  sigma^2 estimated as  0.9222         2. ARMA model              > autofit(t)         $phi         [1] 0              $theta         [1] 0              $sigma2         [1] 0.9130698              $aicc         [1] 279.4807              $se.phi         [1] 0              $se.theta         [1] 0        3. ARIMA model                  > auto.arima(t)         Series: t          ARIMA(0,0,0) with non-zero mean               Coefficients:               intercept                  1.2623         s.e.     0.0951              sigma^2 estimated as 0.9131:  log likelihood=-138.72         AIC=281.44   AICc=281.56   BIC=286.67      The auto.arima function automatically differences time series: we don't have to worry about transformation.              > auto.arima(AirPassengers)     Series: AirPassengers      ARIMA(0,1,1)(0,1,0)[12]                              Coefficients:               ma1           -0.3184     s.e.   0.0877          sigma^2 estimated as 137.3:  log likelihood=-508.32     AIC=1020.64   AICc=1020.73   BIC=1026.39`      Which model should I select to get p,q values & for forecasting purpose?


title: What kind of learning algorithm is suitable for classifying unknown number of groups of unlabeled data?

text: I'm looking for a learning algorithm to analyze data, then group it and define how many classes are fit to distinguish the data set. Type of experiment is in the context of evolutionary biology, how a set of bacteria evolve over the course of time.   1. The data is initially unlabeled so I'm assuming I should look into **unsupervised learning** algorithms.   2. I'm trying to classify instances of the training set into groups, so I thought **clustering** types of algorithms would be more suitable (following a pre-processing of **PCA** for an easier grasp of the data).   3. The number of groups are unknown, so something like K-Means might be more troubling since the number of **centroids** is preselected before the analysis step. I did some reading/research and stumbled upon **Hierarchical Clustering** , which seemed like a great candidate for this problem.


title: Conditional logistic regression model does not converge but logistic regression model does

text: I am running an analysis where I have 2500 cases and 2500 controls. The cases have disease A, and the controls do not. I am trying to see if having disease A increases the odds of various diseases. For the sake of simplicity, we can focus on one disease, call it disease B. D = 1 if disease B present, 0 otherwise E = 1 if disease A present, 0 otherwise I am also including in the model a measure of healthcare utilization. F is a positive integer proportional to an individual's utilization of healthcare. I am running the logistic regression model as such in R:               glm(D ~ E + F, family = "binomial")       Now, this works fine. However, when I try to run conditional logistic regression, it gives me an error:               library(survival)     clogit(D ~ E + F, strata(matched.pairs))     Error in fitter(X, Y, strats, offset, init, control, weights = weights,  :       NA/NaN/Inf in foreign function call (arg 5)     In addition: Warning message:     In fitter(X, Y, strats, offset, init, control, weights = weights,  :       Ran out of iterations and did not converge      I have tried different strata, including dividing the individuals into quantile bins based on F. It does not seem to change anything. (note: pairs are matched on age, gender, race, and F) This occurs only when I run it on a larger sample size. I ran this same analysis on a sample size of 200 (100 cases and 100 controls) and it worked fine. When I use a sample size of 5000, I get the above error. I also made sure that at least 10 cases and 10 controls had the disease in question (disease B, for this example). I am not sure why logistic regression runs fine when conditional logistic regression does not. Can anyone offer me any advice? Thank you all in advance.


title: What is the right attitude toward open source machine learning toolkits?

text: There are lots of machine learning toolkits nowadays, such as _weka_ , _sklearn_ , `R` libs. If we choose to use these toolkits, besides that it is convenient, sometimes we might lose control of what is really happening. For example, in some learning methods the value of features should be scaled to $[0,1]$, but there may be lots of different methods which can achieve this. When we use a scaling function provided by a specific toolkit, we just don't know what's really happening, and this will cause trouble, say, when we want to check if different scaling would affect the performance of learning algorithms. If the learning system is entirely implemented by ourselves, there could not be that much confusion, since we know exactly what happens. So what is the right attitude toward open source machine learning toolkits? Just use them for the first phase, and then implement manually without the help of learning toolkits in order to have more control over algorithms, or can we just rely on such kinds of learning toolkits?


title: Link Functions where MLE=OLS

text: This is a follow-up of a question that I posted previously. I'm trying to get parameter estimates from two different SAS functions (`Proc REG` and `Proc GENMOD`) to equal each other when run on the same data set. `Proc REG` uses OLS/WLS to estimate parameters; `Proc GENMOD` uses MLE. My understanding is that OLS is equal to MLE when the assumed distribution is normal. In my program, I can get the parameter estimates to be equal if I select a normal distribution for `Proc GENMOD` _and if my link function is set to identity_. This works whether I choose to use weights or not. However, if I change my link function in `GENMOD` to `Log`, and change the response variable in `Proc REG` from Y to `LOG(Y)`, the parameter estimates no longer match. Is there something about the use of the log link function that no longer makes OLS=MLE? I've only ever read that equivalence is simply based on the assumption of a normal distribution.


title: Panel data visualization

text: Does anyone know how to graph changes in a dependent variable _over time_ , based on changes in _several_ independent variables (i.e., graphs for panel data regressions)?


title: Good reference book for statistics in medicine

text: I am looking for a good reference book for applying statistics to medical problems. I am interested in a practical book that describes concepts including statistical power, tests for comparing groups, and regression analysis in the context of medical problems.


title: Any technique other than Logistic

text: All my predictors are binary in nature. So far I have been building the model with the logistic fuction. Could anyone suggest any appropriate statistical technique keeping the following points in mind:   1. Dependent variable - binary   2. All Independent variables - binary    3. Response rate - On the lower side (I only include those which have a respose rate higher than 1%) Thanks for your time in advance Kind regards


title: How does PCA improve the accuracy of a predictive model?

text: I've seen in a kaggle challenge about digit recognition someone who used PCA before decision tree or other techniques. I thought it was just for compressing data but he aimed to improve his score. How can PCA improve score in this case ? Is it because there is less overfitting ?


title: How can I compare the effect of an pre-post intervention between two groups with multiple numeric dependent variables?

text: I have collected data from two groups (age 25-30, age 30+). In each group I have 10 patients. I performed several medical tests (i.e., R.B.C, Liver test, etc.) on the patients, before and after treatment. The data look like this:               Group 1                         Test1           ...          Test 10                       Before    After                  Before    After          patient 1       10        16                     13       patient 2       12        11                     19      It's the same for Group 2. ### Question How can I test the effect of the treatment in both groups?


title: How can I estimate a "confidence" in a score / rating

text: I'm sure there's a more technical term for what I'm looking for, but I don't know what it is (it's not confidence intervals, I think) I'm trying to make a system where people get ratings (from 1 to 5 stars) over time, and I want to combine all those ratings into one individual score for each person. Now, this score alone is not a very good thing, I'd like to have a second value that shows how "reliable" / "accurate" this score is. Let's say that:   * If someone gets 1 rating of "3", his score is "3", but we're not very sure of that.   * If someone gets 10 ratings of "3", we're _very_ confident he's a "3". Let's say that more than 10 ratings that are identical are not a stronger signal than 10 ratings, there's a point where we're "confident enough". Let's call this 100% confident.   * If someone gets 2 x"1", and 2 x "5", the average is also a 3, but the person is all over the place, and I want to reflect this by giving this score a _very_ low confidence value. So, I'm pretty much looking for something that reflects variance, but also reflects how many "consistent" ratings someone has got. Or how "stable" he is at his score. Can anyone think of either an algorithm for doing something like this, or point me in the direction of tools that you'd use for such a thing? Thank you!   Daniel


title: Generalized linear models (covariates and splitting files)

text: I am using generalized linear models with one response variable and 6 predictors (1 covariate and 5 factors). I want to assess the effect of smoking on my response variable. When I split my participants in sex groups and then in three smoking groups (smokers, ex-smokers, non-smokers) I find weaker associations than when I add smoking in the model (as a factor since it is a categorical variable) and split the participants in sex groups. Could anyone please explain why is there this difference and what would be the best practice in this model? Thank you! ps: when I split in three smoking groups, the associations are similar per smoking group, they are just weaker.


title: Binning answers from a Likert-scale question

text: I have six questions that each have a five-point Likert scale. We are interested in whether participants, overall, "agree" or "disagree" (or are neutral) toward the questions. To do this, we currently calculate the average of the items, which results in an average from 1 to 5, then we collapse it into 3 bins, where 5 and 4 are "Agree", 3 is Neutral, and 2 and 1 are "Disagree". We then compare this using a $\chi^2$ Goodness of Fit to the expected outcome, which we assume is randomly uniform: we anticipated that as 40% = agree, 20% = neutral, and 40% = disagree, since on a five-point scale 2 options are agree, 1 is neutral, and 2 options are disagree. However I have doubts that we can average on a five point scale, then collapse into a 3 point scale and claim the expected averages are the same. I'm wondering if we can make the expected values claim if we collapse the categories into a 3-point scale before calculating the average. My colleagues don't see the problem. Is there a mathematical argument that shows whether these approaches (collapse before average, or collapse after average) are equivalent?


title: ProjectTemplate equivalent for Python?

text: Is there ProjectTemplate equivalent in python? I.e. a template that defines and helps carry out an efficient workflow for a statistical project? If there is nothing specifically designed for statistics, has anyone successfully used a more generic python project workflow management tool for a statistics project? If so, please give some description of the process.


title: Power estimates for an interaction in a 2X2 factorial design with a proportion outcome

text: I am trying to calculate the power that I will have to detect an interaction between 2 different interventions in a 2X2 factorial design. We aim to have about 850 participants in each of the four boxes, or roughly 1700 per arm for a total of 3400. The study outcome is prevalence at 2 years, with expected prevalence of 2-8% depending on the arm. I was wondering if someone had a good equation to calculate what power we would have to detect an interaction, and what size of interaction? I have not had any luck with an internet search, although have been playing around with factorialsim in stata 12.0, using my proportions (or prevalence) as means, which feels wrong.


title: Plotting Scatterplot Matrix or Correlation matrix or both?

text: I have a problem where I want to use a classifier for it. So I defined a set of features and created a dataset. Now I want to generate some plots to understand the features. I came across the Scatterplot Matrix or the Correlation Matrix, but I'm so confused about if they are the same thing or not. 1 - Visualizing correlation matrices 2- Scatterplot Matrices So my questions:   * Are they the same thing?    * Does plotting both make sense? i.e. each plot convey a different message?


title: Kalman filter transition matrix

text: Hi guys I am trying to writ e a code on python to correct forecast data using Kalman Filter. I am following the equations and recommendations in this link : http://www.swarthmore.edu/NatSci/echeeve1/Ref/Kalman/ScalarKalman.html. At some point, I have written a function which satisfies me since I don't know whatsoever how to use the pykalman module to correct estimations. I am using just one parameter (unidimensional and univariate time series), so my matrices form will be a numerical value. My main problem is what is the transition matrix? In many paper they just use a transition matrix of 1 but I think this is not rigorous since it would mean that: forecast(t) = A * forecast(t-1) We do not have a constant value of forecast so what would be my transition matrix? Hypothesis 1: Can I just try to see a correlation coefficient between two estimated values then take the mean of those coefficients for my first day of forecast.By the way this transition matrix will be just an initialization and will change for each forecast following


title: When is it appropriate to use PCA as a preprocessing step?

text: I understand that PCA is used for dimensionality reduction to be able to plot datasets in 2D or 3D. But I have also seen people applying PCA as a preprocessing step in classification scenarios where they apply PCA to reduce the number of features, then they use some Principal Components (the eigenvectors of the covariance matrix) as the new features. My questions:   * What effects does that do to the classification performance?   * When to apply such a preprocessing step?   * I have a dataset with 10 features as real numbers and 600 binary features that represent categorical features, using one-to-many encoding to represent them. Would applying PCA here make sense and make a better results? p.s. if the question is too broad, I would be thankful if you provide a paper or tutorials that explains better the details of using PCA in that manner. p.s. after reading a little, i found that it could be better to use Latent Semantic Analysis to reduce the number of binary features for the categorical attributes? So I don't touch the real-valued features, but only preprocess the binary features and then combine the real-valued features with the new features and train my classifier. What do you think?


title: Does the sequence satisfy WLLN?

text: Could you help me prove that the following sequence of independent random variables satisfy Weak Law of Large Numbers? $$P(X_n = n - \sqrt{3}) = P(X_n=-n- \sqrt{3}) = \frac{1}{2} $$


title: How to use statistics to learn about voting?

text: I have a data set with 1500 individuals who were surveyed about many things including demographics, attitudes and opinions about politics and economic status. The most important question in the survey was about vote intention on the next election for president. The question had the four candidate options, plus options for 'I don't know yet', and for 'I'm not going to vote'. Among the 4 candidates, 2 are strong, and 2 are very weak. The question I'm trying to answer is: among the undecided (the ones who say I don't know yet), what are the characteristics of the voters who are more likely to vote for each of the candidates? For example, are male or female among the undecided more likely to vote for candidate A? How can I do that in R?


title: Interpretation of PCA biplot?

text: I just ran my first ever PCA, so please excuse any naivety on my part. As input, I used five years worth of the following:   * S&P;/ASX 200 A-REIT   * S&P;/ASX 200 Consumer Discretionary   * S&P;/ASX 200 Consumer Staples   * S&P;/ASX 200 Energy   * S&P;/ASX 200 Financial-x-A-REIT   * S&P;/ASX 200 Health Care   * S&P;/ASX 200 Industrials   * S&P;/ASX 200 Information Technology   * S&P;/ASX 200 Materials   * S&P;/ASX 200 Resources   * S&P;/ASX 200 Telecommunication Services   * S&P;/ASX 200 Utilities Using R, I simply ran the following commands:               arc.pca1 <- princomp(sp_sector_data, scores=TRUE, cor=TRUE)     summary(arc.pca1)     plot(arc.pca1)     biplot(arc.pca1)      ### Summary               Importance of components:                              Comp.1     Comp.2     Comp.3     Standard deviation     2.603067 1.05203261 0.88394057     Proportion of Variance 0.564663 0.09223105 0.06511258     Cumulative Proportion  0.564663 0.65689405 0.72200662                                     Comp.4     Comp.5     Comp.6     Standard deviation     0.84122312 0.76978259 0.73901015     Proportion of Variance 0.05897136 0.04938044 0.04551133     Cumulative Proportion  0.78097798 0.83035842 0.87586975                                     Comp.7     Comp.8     Comp.9     Standard deviation     0.66409102 0.62338449 0.52003850     Proportion of Variance 0.03675141 0.03238402 0.02253667     Cumulative Proportion  0.91262116 0.94500518 0.96754185                                    Comp.10    Comp.11      Comp.12     Standard deviation     0.45637805 0.42371864 0.0409804189     Proportion of Variance 0.01735674 0.01496146 0.0001399496     Cumulative Proportion  0.98489859 0.99986005 1.0000000000      ### Loadings               Loadings:             Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7     RE      -0.235         0.520 -0.533 -0.438  0.355 -0.150     disc    -0.332               -0.125                0.294     staples -0.295  0.226                      -0.211  0.554     energy  -0.332 -0.251         0.172  0.176        -0.130     fin_RE  -0.323               -0.118 -0.130         0.384     health  -0.224  0.465 -0.124 -0.193  0.603  0.537 -0.112     ind     -0.337                                               IT      -0.224  0.145 -0.757        -0.461        -0.312     mat     -0.329 -0.351         0.295         0.126 -0.116     res     -0.335 -0.350         0.297         0.123 -0.133     telco   -0.161  0.609  0.327  0.609 -0.311        -0.113     util    -0.270  0.160  0.146 -0.256  0.234 -0.694 -0.509                  Comp.8 Comp.9 Comp.10 Comp.11 Comp.12     RE      -0.217                                    disc     0.309  0.567  0.596                      staples -0.688        -0.141                      energy         -0.215  0.240  -0.783  -0.165      fin_RE   0.374 -0.724          0.207              health                                            ind      0.398  0.311 -0.743  -0.221              IT      -0.183                                    mat     -0.127                 0.461  -0.638      res     -0.123                 0.226   0.752      telco    0.116                                    util                           0.115                                  Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6     SS loadings     1.000  1.000  1.000  1.000  1.000  1.000     Proportion Var  0.083  0.083  0.083  0.083  0.083  0.083     Cumulative Var  0.083  0.167  0.250  0.333  0.417  0.500                         Comp.7 Comp.8 Comp.9 Comp.10 Comp.11     SS loadings     1.000  1.000  1.000   1.000   1.000     Proportion Var  0.083  0.083  0.083   0.083   0.083     Cumulative Var  0.583  0.667  0.750   0.833   0.917                         Comp.12     SS loadings      1.000     Proportion Var   0.083      ### Scree Plot ![enter image description here](http://i.stack.imgur.com/6edCJ.png) ### Biplot ![enter image description here](http://i.stack.imgur.com/QrCyS.png) ### Is this useful? Am I right in assuming that these indices are correlated with each other? Does the biplot show some sort of clustering? What if anything, does any of this mean?


title: What can I read to give me a meta-view of statistics as a field

text: What can I read that will give me a meta-view of the diverse field of statistics and data science? With few exceptions much of what I get my hands of goes straight into formulae and methodologies. Preferably something sufficiently high level as to bring in diverse areas such as econometrics, psychometrics, machine learning, etc. To be more concrete, preferably something that:   * Summarises and discusses the various branches of statistics/data science; what problems are encountered in each branch.   * Talks about differences between the branches and their histories.   * Contrasts methodological approach.


title: expectation of multivariate discrete distribution

text: The expectation of a function $f(x)$ over a probability distribution $p(x)$ where $x=(x_1,\dots,x_n)$ and $x_i \in \left\{1,\dots,K\right\}$ requires a summation over all possible $K^n$ combinations. $ E[f(x)] = \sum_x f(x)p(x)$ For large scale systems this exponential grow up is formidable. I am wondering whether it is possible to approximate this expectation somehow. For a simpler case, which I am mainly interested in; consider a fully factorized distribution $p(x) = \prod_{i=1}^n p(x_i)$. For this case is there any efficient way of calculating the expectation? $ E[f(x)] = \sum_x f(x) \prod_{i=1}^n p(x_i)$ If $f(x)$ was also $\prod f_i(x_i)$ than things would be simple due to $E[XY]=E[X]E[Y]$ however for general nonlinear $f(\cdot)$ things get coupled. For the simplest case which is still of interest to me consider multivariate Bernoulli variables. Again assume $p(x) = \prod p(x_i)$ and $p(x_i=1)=\alpha_i$. Is there a cheap way of calculating $ E[f(x)] = \sum_{x \in \left\{0,1\right\}^n} f(x) \prod_{i=1}^n p(x_i)$ As I mentioned, the result does not need to be exact, a fast approximation scheme is what I am looking for. The first approximation idea that comes to my mind is Monte-Carlo approximation; however I don't think that it will be efficient and accurate enough.


title: Notational question about Bayesian prior

text: This is a simple question about notation. Let $\mathbf{p}$ be a Bayesian prior. Let $\mathbf{x}$ be the data, which is a random vector. Let $\mathbf{S}$ be some indexing set for $\mathbf{x}$. Then I've used the notation $\mathbf{x}_{\mathbf{S}}$ for the restriction of $\mathbf{x}$ to $\mathbf{S}$, and $\mathbf{p}_{\mathbf{S}}$ for the marginal prior of $\mathbf{p}$ on $\mathbf{S}$. I've been using $\mathbf{y}$ as shorthand for $\mathbf{x}_{\mathbf{S}}$. What I was wondering was, what would be suitable shorthand for $\mathbf{p}_{\mathbf{S}}$. Is it permissible (and not abuse of notation) to just use $\mathbf{p}$, or should I use a different notation like $\mathbf{q}$? It seems to me that using the full prior $\mathbf{p}$ in the context of the subvector $\mathbf{y}$ should not cause any confusion, but I don't want to use incorrect notation.


title: How do I estimate a differences in differences model when the dependent variable has many zeros?

text: Is there any way to run an OLS difference in differences model when the dependent variable (investment) has lots of observations which are truly zero? I don´t know how to add clarifications. My problem is indeed the normality assumption. I am using firm level panel data. Many firms (around 60%) do not invest each year. That`s why I have so many zeroes.


title: Model some logarithmic looking data in R

text: I have the following data in R               x <- c(0.1,0.2,0.3,0.6,0.8,0.9,1)     y <- c(90,96,97.7,99.3,99.65,99.95,100)      I'm trying to find a logarithmic equation that best fits these points. I'm not sure what the equation would look like, but probably something like one of these   * a*log(x)+b   * (a*log(x)+b)/(log(x)+c)   * a*log(b*x)+c   * etc What kind of curve do you think best fits this data. And how can I find out?


title: What is the maximum value in a finite selection of a normally distributed variable?

text: A parameter of an object is normally distributed with a mean `m` and a std. dev. `s`. If `r` such objects are randomly selected, what is the maximum **expected** value `M` of the parameter in the selection? Cross-posted from: http://mathoverflow.net/questions/105903/what-is-the- maximum-expected-value-in-a-finite-selection-an-object-with-normal-di


title: Distance invariant approaches to find the "main" difference between two distributions in R$^n$

text: If I have a population of vectors in R$^n$ and some special subset has a different distribution, I could try and use PCA to describe the main axes of these distributions and, if they are aligned differently, that would give me a quick way to describe the difference in orientation between the two. However PCA is dependent on the scale of the components of the vectors. So what is an analagous approach that can cleanly describe the "main" difference in the the two populations given that they have different covariance matrices, but doesn't depend on the scaling given to the componenents? It is possible the answer is more easily stated in terms of the two Mahalanobis distances one gets from the two distributions. Anyway, I'm sure there is a standard way to describe the "main" difference between two such distributions, but I don't know what it is. P.S. Yes, assume they have the same mean.


title: What statistical tools would you use for identifiying important factors/variables in an experiment?

text: I conducted a experiment where I measured bacteria on a 1000 peoples' hands after performing a repeatable task. ![enter image description here](http://i.stack.imgur.com/8wbDy.png) This scatter graph shows 2 experiments where the humidity of the hands is changed. I then repeated the experiment under 4 different conditions. What would you use to work out whether a variable was statistically important in alterring the results? EDIT: The experiment is actually numerical in nature: The black points are produced by a function: $$Y=f(A,V,\lambda)=\sum_{i=1}^nA_iV_i\lambda_i$$ Where $A$ $V$ and $\lambda$ are drawn from empirical distributions the red points conform to: $$Y=g(A,V,\lambda,\beta)=\sum_{j=0}^i \left(\prod_{k=j+1}^i \beta_k\sigma_j\right)$$ where $\sigma=A\lambda V$ and $\beta$ is an empirical distribution


title: Fixed effects interpretation

text: If I have a fixed effects model where data is separated by states over time, how would the coefficient interpretation on variables differ from an OLS interpretation? And how would I interpret the intercept?Thank you.


title: Can I run a GLMM model when I have one observation for most subjects?

text: I have a binary DV and my panel data set contains more than one observation for only 20% of the subjects which makes it very unbalanced. Is there anything methodologically wrong with doing a mixed model with data that looks like this,i.e. with many subjects having only one observation? I'd appreciate any hints. > No. of observations per subject:     1       2         3     4     5   >  Frequency:                                 4145    821    176    57    19 I am reading "Linear Mixed Models A Practical Guide Using Statistical Software" by West, Welch and Galecki, and my data set does not seem to exactly fit into their description of three types of data they discuss for mixed models, namely Clustered Data, Repeated-Measures and Longitudinal Data.


title: What common forecasting models can be seen as special cases of ARIMA models?

text: This morning I woke up wondering (this could be due to the fact that last night I didn't get much sleep): since cross-validation seems to be the cornerstone of proper time-series forecasting, what are the models I should "normally" cross-validate against? I came up with a few (easy) ones, but I soon realized they were all but special cases of ARIMA models. So I'm now wondering, and this is the actual question, what forecasting models does the Box-Jenknins approach already incorporate? Let me put it this way:   1. Mean = ARIMA(0,0,0) with constant   2. Naive = ARIMA(0,1,0)   3. Drift = ARIMA(0,1,0) with constant   4. Simple Exponential Smoothing = ARIMA(0,1,1)   5. Holt's Exponential Smoothing = ARIMA(0,2,2)   6. Damped Holt's = ARIMA(0,1,2)   7. Additive Holt-Winters: SARIMA(0,1,m+1)(0,1,0)m What else can be added to the previous list? Is there a way to do moving average or least squares regression "the ARIMA way"? Also how do other simple models (say ARIMA(0,0,1), ARIMA(1,0,0), ARIMA(1,1,1), ARIMA(1,0,1), etc.) translate? Please note that, at least for starters, I'm not interested in what ARIMA models **cannot** do. Right now I only want to focus on what they **can** do. I know that understanding what each "building block" in an ARIMA model does should answer all of the above questions, but for some reason I have difficulties figuring that out. So I dediced to try a "reverse engineering" kind of approach. Any input is appreciated, I hope I stated my question clearly enough.


title: How to forecast based on aggregated data over irregular intervals?

text: I am trying to forecast the sales of products in vending machine. The problem is that the machine is filled at irregular intervals and at every fill we only can record the aggregated sales since the last fill of the machine (i.e. we don't have daily sales data). So basically we have data for aggregated sales at irregular intervals. The intervals usually are between 2 days and 3 weeks. Here is example data for one vending machine and one product:               27/02/2012 48     17/02/2012 24     09/02/2012 16     02/02/2012 7     25/01/2012 12     16/01/2012 16     05/01/2012 16     23/12/2011 4     16/12/2011 14     09/12/2011 4     02/12/2011 2      Our current naive algorithm is to calculate average sales per day by dividing the total quantity sold during the last 90 days by 90. Do you have any idea how to improve the forecast of the sales per day? I need to forecast what will be sold at the next visit of the machine. Is it possible to use some sort of exponential smoothing algorithm given the nature of our data? Thanks in advance! UPDATE: Thanks a lot for all the answers and comments. Let me try to give a bit more context (the business case behind the question - very simplified of course). We have hundreds of vending machines. Every day we have to decide which 20 of them to visit for refill. In order to do so we are trying to predict what is the current status of the machines and to select the "emptiest" 20 machines. For each machine and product we are calculating the average sales per day (SPD) using the naive algorithm described above. Then we multiply the SPD by the number of days since the last fill of the machine and the result is the predicted quantity sold.


title: Multiple curves when plotting a random forest

text: I performed a random forest using the `randomForest` package. I know that if I plot the random forest using the `plot()` command, I should get back a graph with number of trees on the x-axis, and estimated error rate on the y-axis, but when I actually plot it, I get multiple curves, all of them start quite high, as expected and drop as the number of trees increase. Well, what do they mean? I assume they are just different measurements, e.g. one is the `oob error rate`, another is the `misclass error rate` and so forth, but how can I extract that information and/or filter the undesired curves?


title: Use of PCA analysis to select variables for a regression analysis

text: I have too many environmental variables to use in a multiple regression analysis. If I use all the variables the models are just too complex. The use of the PCA axes in the regression analysis was impossible to interpret (since there wasn't a clear correlation with environmental variables), so we chose to select a limited number of variables, namely those that had the higher explanation in PCA. A PCA was used for each set of environmental variables, namely the variables related with the structure of the stream, the evolving vegetation, the climatic, the physical-chemical characteristics of the water from summer period and from winter period, separately. The PCA was performed using the correlation matrix option, using the software PC-ORD, v. 4.21 (McCune & Mefford 1999). For each set of variables, only the variables with coordinates higher than 0.20, for the two first axes, of the PCA, were selected to be used in multiple regression analysis. I could not find literature that confirms that it's OK to do this, but I think it is not wrong.


title: Understanding the effect of a continuous random factor in a mixed effects model

text: I understand the effect of a categorical random effect on a mixed effects model in that it performs a partial pooling of the observations by level in the random effect, effectively assuming that the observations aren't independent themselves but only their partial pools are. Also to my understanding, in such a model observations sharing the same random effect level but differing in their fixed effect level will outweigh observations differing in both their random effect and fixed effect levels. What is the effect of a continuous random factor then? Given that a model without the random effect showed that the fixed effect had an effect size X. Should I expect that if the observations in the different levels of the fixed effect came from far ends of the random effect continuum the effect size will become smaller in a model which included the random factor, while if observations in different fixed factor levels had similar random effect values then the effect size would increase?


title: How to fit (linear regression) a ratio of all independent variables?

text: I am trying to find the best fit between an species dataset and prevailing climatic conditions, in order to be able to predict the environmental conditions from the species dataset (paleoclimate research). I have 15 species(sp1-sp15), expressed as relative amounts (some are 0). I have done some data exploration in excel, and have seen that I get good fits using the ratio sp2/(sp1+sp2). I have done multiple linear regression on this dataset, but none give a correlation as good as the ratio I stumbled on. The ultimate goal is of course to test whether this ratio I stumbled upon is yields the best model to explain changes in the climate variable. I would thus like to create the ratios: `sum(spj)/sum(spj)`, where spj refers to one of the species. Here, both the numerator and denominator can be any possible combination of all species. Can I generate a formula model to include a ratio? Can I then select the best model using regbsubsets? **EDIT** Based on answers below, this will not be possible for all 15 variables. How would it be possibel for 5 variables (I have seen this published before, using R). **EDIT** On stackoverflow (where I posted the same question), I got the comment to use the function (I), to write expressions inside formulae. However, if I use               leapsMAT<-regsubsets(x1 ~ I(1+ (y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8 + y9 + y10 + y11 + y12+ y13 + y14 +y15 )/(y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8 + y9 + y10 + y11 + y12+ y13 + y14 + y15)), force.in= FTFALSE,nvmax=15,data=my.data, nbest=1)      the regsubsets doesn't recognize the variables as variables, of course. How do I implement the I(function), or is there another way to deal with this? At the moment I have obtained the same linear regression using two methods:               library(leaps)     attach(my.data)     FTFALSE<-c(FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE)     leapsMAT<-regsubsets(x1 ~ 1+ y1 + y2 + y3 + y4 + y5 + y6 + y7 + y8 + y9 + y10 + y11 + y12+ y13 + y14 + y15 , force.in= FTFALSE,nvmax=15,data=my.data, nbest=1)      And, using a longer method specified here (implemented from somewhere else on the internet):               allModelsList <- apply(regMat, 1, function(x) as.formula (paste(c("x1 ~ 1", namevar2[x]),collapse=" + " )))     allModelsList     warnings()          #Calculating the model     my.data     allModelsResults <- lapply(allModelsList, function(x) lm(x, data=my.data))     allModelsResults          dfCoefNum   <- ldply(allModelsResults, function(x) as.data.frame(t(coef(x))))          dfStdErrors <- ldply(allModelsResults, function(x) as.data.frame(t(coef(summary(x))[, "Std. Error"])))               dftValues   <- ldply(allModelsResults, function(x) as.data.frame(t(coef(summary(x))[, "t value"])))          dfpValues   <- ldply(allModelsResults, function(x) as.data.frame(t(coef(summary(x))[, "Pr(>|t|)"])))      dfpValues     (warnings)     # rename DFs so we know what the column contains     names(dfStdErrors) <- paste("se", names(dfStdErrors), sep=".")     names(dftValues) <- paste("t", names(dftValues), sep=".")     names(dfpValues) <- paste("p", names(dfpValues), sep=".")          # p-value for overall model fit     calcPval <- function(x){         fstat <- summary(x)$fstatistic         pVal <- pf(fstat[1], fstat[2], fstat[3], lower.tail = FALSE)         return(pVal)     }          # Before creating ONE data frame with all important entries,     # we need to compute some more indices      NoOfCoef <- unlist(apply(regMat, 1, sum))     R2       <- unlist(lapply(allModelsResults, function(x)                               summary(x)$r.squared))         adjR2    <- unlist(lapply(allModelsResults, function(x)                                   summary(x)$adj.r.squared))     RMSE     <- unlist(lapply(allModelsResults, function(x)                               summary(x)$sigma))     fstats   <- unlist(lapply(allModelsResults, calcPval))                    # now we can combine all the data into one data frame     results <- data.frame( model = as.character(allModelsList),                            NoOfCoef = NoOfCoef,                            dfCoefNum,                            dfStdErrors,                            dftValues,                            dfpValues,                            R2 = R2,                            adjR2 = adjR2,                            RMSE = RMSE,                            pF = fstats  )     results[1:20,]     # round the results     results[,-c(1,2)] <- round(results[,-c(1,2)], 3)     results          model.maxRadj<-which(results$adjR2 == max(results$adjR2), arr.ind = TRUE)     maxRadj<-results[model.maxRadj,]      Many thanks in advance! Please let me know if more information is required.


title: Using R for an artificial neural network (ANN)? Or other frameworks?

text: I'm planning on starting a ANN classification project pretty soon. It will use 15 input variables (medical data with 25k records) to predict an outcome. I have been following the ML course on Coursera, where Octave is used, but I'm used to `R`. I've been reading about ANN in `R` and I got the feeling the solutions are not optimal. Main source of worry: this thread. Here it is said that `R` implementations of ANN are not great, and other frameworks are 'better'/faster. That's why I've been looking into Encog, and Encog Workbench, which is pretty neat. But still, using R would be nice. I could also consider using Octave. Can you comment on the quality of the `R` ANN frameworks? Would you advice me to still use `R`? I really need the best accuracy I can get, as these data are most likely not be perfect.


title: Conceptual question on pattern recognition

text: I am just a beginner in machine learning and finding tough time to understand few concepts. I have 10 strings of characters of length 50 each. These are the training data $T(i)$ for $i=1:10$. The first 5 come from class $A$ and next 5 from class $B$. Then there is an incoming string\test Test of same length. I calculate the number of common characters between each $T(I)$ and Test. Let this be a distance measure $dist(T(i),Test)$. Test will belong to or match a $T(i)$ if the number is maximum i.e Test $\epsilon$ max $dist(T(i),Test)$   * Can this distance count serve as a feature for recognizing task maybe feeding this distance into a Bayes classifier or k-nn?   * How can I apply clustering or classification for this kind such that I may know the incoming string is most similar to the training string and assign it to that class? I know it is vague or maybe inappropriate but can somebody explain with a code how to do this and if it is valid to do so?


title: Proof of how the normal distribution arises as the limit of the binomial — where do $\pi$ and $e$ come from?

text: I've always thought the emergence of the normal distribution was kind of magic, specifically that $\pi$ and $e$ both emerge in the distribution, even though they don't exist in the binomial. I've googled heavily, but haven't found a very good guide showing how these constants arise, but would like to see it. I'm horrible at proofs, so can someone walk me through taking binomial to the limit to arrive at the normal? Wolfram says de Moivre developed this before 1783. Also, if possible, please don't use the CLT for the proof. Perhaps I'm wrong, but my understanding is that the CLT wasn't proven until the 19th century.


title: Feed forward neural networks for the analysis of censored survival data, how should I implement it in R?

text: Recently, I have read an article which name is “ _Feed forward neural networks for the analysis of censored survival data: A partial logistic regression approach_”. Without a math background, I catch only a little about the ANN applying the analysis of censored survival data. In the training group, each subject is replicated for all the intervals in which the subjects is observed and coupled with the event; why is this done? Especially, in the testing group, indicator each subjects are replicated into full number of time interval of observed with all event indicator as zero? I know this process can take into account the censoring data. In addition, are the output of neural networks are conditional probabilities of failure? My problem is how to produce the survival curve using the output data. Does anybody know the R code or MATLAB code to perform this whole process? Or give me some suggestion to find answers! The following R code is my try on this method but I can't go on it, for I don't know draw the survival curve depending on the output conditional failure probabilities!               dat<-read.csv("traininglj.csv",header=T)     tt <- dat$time <- as.numeric(cut(dat$TTR,c(0,6,12,18,24,30,36,42,48,54,200)))     dat2 <- dat[rep(1:nrow(dat), tt), ]     time2 <- NULL     for (i in 1:length(tt)) time2 <- c(time2, 1:tt[i])     dat2$time <- time2         dat2$Recurrence <- 0     dat2$Recurrence[cumsum(tt)] <- dat$Recurrence     write.csv(dat2,file="result.csv")     mydat <- apply(dat2[,13:23],2,function(x)(x-min(x,na.rm=TRUE))/         (range(x,na.rm=TRUE)[2]-range(x,na.rm=TRUE)[1]))     training <- cbind(dat2[,1:12],mydat,dat2[24])     library(nnet)     library(lattice)     attach(training)     dat.net <- nnet(Recurrence ~ time+ALT+ALB+PLT+INR+age+MELD+logAFP+Diameter             +sex+number+Gstage+HBsAg+AN+MVI,          data = training,          size = 12,          decay=0.025,         maxit = 1000,         entropy=TRUE,         trace=TRUE)


title: How to enhance the performance of naive CMDS?

text: I'm performing Classical MDS on the distance over columns of binary matrix. The result is like this: ![enter image description here](http://i.stack.imgur.com/XkB9J.png) The points lies on two lines vertical to each other. I don't think useful information is shown in this MDS chart. So I'm looking forward to the ways to enhance the classical MDS perfromance. The data source is a 100 * 600, {0,1} valued matrix. I use normalized euclidean matrix between columns.               dist(col1, col2) = EuclideanDist(col1, col2) / sqrt(Norm2(col1)*Norm2(col2))


title: How to perform cross validation of a model predicting categorical variable in SAS?

text: I have a prediction outcome that is dichotomous while all other predictor variables are binary and categorical. How can I perform cross validation in SAS and see the ROC curve? I need to create bayesian network.


title: Web World Wide Data Mining

text: I have been working on a project out of curiosity to get millions of URLs returned by search engines and save them to a text file or any other data source of mine. with time my I started my research by looking at Google Search API old one and new ones but they just returns you couple of hundred of URLs only, I also tried looking at different other search engines but couldn't find much. Now since I am more curious I started looking at how search engine works, how they index million of websites, "web crawler" was the answer but am still looking for the best way to get URLs of websites for a specific search query using web crawler might take me ages to get URLs I want. Since this question is related to Data Mining so posting here, Would really appreciate if someone guide me in right direction please. **what I actual need is "Millions of URLs for a specific query" e.g. all websites returned for "Apple"** I don't mind reading books at all, and also if you can state that working on such kind of projects can help me in any type of career. Thanks for your time


title: Self organizing maps vs. kernel k-means

text: For an application, I want to cluster data (potentially high dimensional) and extract probability of belonging to a cluster. I consider at the moment Self organizing maps or kernel k-means to do the job. What are the pros and cons of each classifier for this task? Am-I missing others clustering algorithms that could be performant in this case?


title: Probability of a given password

text: A password consists of 4 alphabet letters and 4 numbers. Calculate the following two probabilities:   1. $p_1$: the probability that the letters are all equal and that the numerical part contains one eight.   2. $p_2$: the probability that the password has 3 numbers followed by 4 letters. Although it sounds like an easy question, but how would I apply the definition of permutations and combinations here? Here is how I thought of solving it. $p_1= (1/21)^4*(1/10)*(9/10)^2$ Do I need to calculate all the possible combinations here? $p_2= 1/\binom{7}{7}=1/7!/(7-7)!= 1/7!$ since we are considering only one case among a permutation of 7 elements over 7 places.


title: Is a mixed/ random effects model required if fixed effects model shows no pattern in the residuals?

text: I've got some data on discrete flood events (response variable - duration) on several rivers. I've modelled the response variable by site and time using a Generalized Linear Model:                mod = glm(duration ~ site + time)      I've then looked at the residual plots and cannot find any evidence for temporal correlation in the residuals. Does this mean that this model approach is valid (and that using a mixed/ random effects model is unnecessary) (or should one always use mixed/ random effects models for time-series data)? Any links to reading references regarding this subject would be brilliant,


title: Are empirical risk minimization and M-estimators the same thing?

text: Would it be true to say that Empirical Risk Minimization and M-Estimation are the same thing?


title: Probability density function from a set of data

text: Given a set of experimental data f={fk} for k=1,2,...,N, how can I find the best pdf which approximate the set fk? In the case of a gaussian distribution, applying the maximum likelihood criteria, it's easy to find the couple (μ,σ), but if I don't know anything about the distribution of the data, is the maximum likelihood estimate still valid? And if the answer is 'yes' what is the formula? Many thanks.


title: Should I interpret based on descriptive or inferential statistics in the present case?

text: In our study which is completed, I have three groups: Drug A (n=87), Drug B (=62) and Drug C (n=5). The Shapiro-Wilk test, histograms and side-by-side box plots were used to assess normality of the data of the individual groups with respect to all the parameters used separately (weight, BMI, blood sugar, total cholesterol, triglycerides, LDL-cholesterol, HDL-cholesterol, clinical / illness / adverse effect / quality of life, etcc, scores). So, ultimately we check for non-normal distribution and outliers, and are concerned about unequal sample sizes (n=5 in the Drug C group), in addition to looking at descriptive statistics (mean ± SD and proportion of % wherever applicable), nonparametric tests were carried out for the inference: For between-group comparisons, Kruskal Wallis H and subsequently for the therapeutic effect where p<0.05 obtained in Kruskal Wallis H, Mann Whitney U test as post-hoc test were used. And in case of within-group comparisons, Freidman’s test (for five evaluable visits data) and subsequently for p<0.05 for the therapeutic effect where p<0.05 was obtained in Freidman’s test, Wilcoxon-signed rank test as post-hoc test were done and p<0.05 was obtained using Wilcoxon-signed rank test in two (out of three) groups. Also, to ascertain the between-group significance difference, Mann Whitney U test was applied. Kindly advise me:   1. First of all, have I made the right choice of tests and am I going in right direction?    2. With regards to the present analytical sample of N=154 in total, for arriving at final conclusion should I depend on descriptive or inferential statistics or both?    3. In the case where the mean ± SD changes and proportion of percentage (%) are found comparatively higher in case of Drug C which has sample size n=5 compared to other two drug groups but in case of between-group and/or within group inferential statistical analysis, (wherever) found non-significant difference / change in case of Drug C, can it be ‘said’ superior to other two drugs ‘on the basis of descriptive statistics’ or can I at least make a statement that mean ± SD and proportion of % changes at endpoint (or at other time-point) compared to baseline were found higher in case of Drug C compared to other two or other one? So, what should be the (cautious) interpretation / statement in case of Drug C compared to other two drugs?


title: Convergence problem

text: I practiced some problem (not homework) as below. Suppose that $X_1,X_2,\cdots$ are independent and identically distributed real-valued random variables with $E|X_1|=\infty$. Show that $\sum_{n=1}^\infty P({|X_n| \geq a\cdot n})=\infty$ for each $a >0$. I obtained that $\infty = E(|X_n|) = \int_{0}^\infty P(|X_n| >x)dx \leq \sum_{n=0}^\infty P(|X_n|>n)$, but I couldn't get the final step $\sum_{n=1}^\infty P({|X_n| \geq an})=\infty$ for each $a >0$ ( or can I use the statment to imply the result?) Please provide some hints, thanks!


title: What is the ANOVA model to use for a 2x2 factorial experiment with units measured 3 times each?

text: If I were to set up an experiment examining the effects of exercise, and coffee consumption on the number of hours slept, what would the be the ANOVA model that I should use, given:   * 16 volunteers randomly assigned to two groups:      * 1st group will drink cup of coffee in the morning,      * 2nd group nothing   * Within each group:       * 4 volunteers will exercise once per day,      * the rest will not.   * Each subject will be tested on n=3 separate occasions, and the number of hours slept will be recorded. I find it extremely confusing to get the correct model.


title: Large data set for testing kernel logistic regression

text: Does anyone know of a large data set (upwards of $10^7$ rows, but I'll take $10^5$ as well) that would be appropriate for testing kernel logistic regression? Continuous variables, 2 to 50 independent variables? Would prefer business like data as aposed to scientific (weather, chemical, bio ...).


title: How to export and use results of PCA from R?

text: My ultimate goal is to run a cluster analysis on a data set with > 1 million records. The input variables for the cluster analysis will be the results of a Principal Component Analysis, as well as other variables not included in the PCA, for a total of maybe 10 variables input into the clustering (the variables I input into the PCA were all very highly correlated with one another while the other variables are not so I chose not to include them in the PCA).               #read data     mydata <- read.csv('mydata.csv')           #import library for robust methods because my data contained outliers     library(rrcov)           #run robust PCA method called PcaCov     pcaR <- PcaCov(~., mydata, na.action=na.omit, center=TRUE, scale = TRUE, k=8)          #look at results     summary(pcaR)     screeplot(pcaR)     pcaR@loadings      From the results, I have decided I would like to retain the first three components, which capture ~87% of the total variance in the dataset. Now I want to extract/save/export these first three components for use in the cluster analysis with my other variables. How do I do this?


title: Why is the R2 or adjusted R2 not used without a constant?

text: The title basically says everthing of importance. I read several times that a R2 is not used/interpreted if one is not using a constant. My question is: why? And if this statement is valid, I would be glad for valid source because I couldnt find one myself.


title: PCA of non-Gaussian data

text: I have a couple of quick questions about PCA:   * Does the PCA _assume_ that the dataset is Gaussian?    * What happens when I apply a PCA to inherently non-linear data? Given a dataset, the process is to first mean-normalize, set the variance to 1, take an SVD, reduce rank, and finally map the dataset into the new reduced- rank space. In the new space, each dimension corresponds to a "direction" of maximal variance.   * But is the correlation of that dataset in the new space always zero, or is that only true for data that is inherently Gaussian? Suppose I have two datasets, "A" and "B", where "A" corresponds to randomly sampled points taken from a Gaussian, while "B" corresponds to points randomly sampled from another distribution (say Poisson).   * How does the PCA(A) compare to the PCA(B)?    * By looking at the points in the new space, how would I determine that the PCA(A) corresponds to points sampled from a Gaussian, while PCA(B) corresponds to points sampled from a Poisson?    * Is the correlation of the points in "A" 0?    * Is the correlation of points in "B" also 0?    * More importantly, am I asking the "right" question?    * Should I look at the correlation, or is there another metric that I should consider?


title: Fitting distributions and comparing tails

text: Suppose I have some data points (0.6695, 0.5968, 0.7641, 0.7252, and 0.7779) and want to fit different distributions to this dataset, say, log-normal and log-t (with different _df_ ), and then compare the tails of these distributions using graphics. How this can be done in R?


title: Which distributions have closed-form solutions for maximum likelihood estimation?

text: Which distributions have closed-form solutions for the maximum likelihood estimates of the parameters from a sample of independent observations?


title: Interpret significant interaction but no main effect in repeated measures ANOVA

text: I have a model in which two groups are measured on one DV at pre- post- and 3 months following an intervention (the control group did not receive the intervention). I have an interaction effect, which I believe implies that the combination of group and time together influence the DV, and yet the same SPSS ANOVA output shows no main effects. How should I interpret the data? Is there an easy way to do further analyses with SPSS?


title: Testing for equality of variance - with repeated measures?

text: I have 40 "subjects", each of them did 5 types of exams. The question I wish to test is if there is a difference in variance between the exams. That is, I would like to run something like a "Bartlett's test", but for "repeated measures". Does such a test exists? Thanks.


title: Given MCMC samples, what are the options for estimating posterior of parameters?

text: Markov chain Monte Carlo (MCMC) is a class of algorithms for sampling from probability distributions. In the end, given a parametric model explaining the data, MCMC could also be used for parameter estimation, in which case I wonder what the different options for processing MCMC samples are. For monomodal distributions for instance, I can think of averaging the $n$ last samples of a run, averaging the last sample of different MCMC runs. What are the options for estimating posterior of parameters?


title: Noisy linear relationship: Can the functional form be known?

text: Let's say I know the relation between `x` and `y` is linear yet noisy. Given a noisy `(x,y)` dataset is there a way to deduce if it was most likely generated out of an underlying function (a) $y = a + b \times x + noise$ or via (b) $x=c + d \times y + noise$ I'm not sure if I know the nature of the noise but let's assume it is gaussian and non correlated. Normally I would use regression to determine the coeficents and I do know that the estimated parameters differ depending on whether `x` is taken to be the independent variable or `y`. **But I don't know if there is a statistical way to differentiate the two cases (a) and (b) that I described.**


title: How can I find $\text{Cov}(X_k,X_5)$

text: > There are $10$ empty boxes numbered $1, 2, \dots , 10$ placed sequentially > on the circumference of a circle.We perform $100$ independent trials. At > each trial, one box is selected with probability $\dfrac{1}{10}$ and one > ball is placed in each of the two neighbouring boxes of the selected one. > > Define $X_k$ to be the number of balls in the $k$-th box at the end of $100$ > trials. > > (a) Find $E[X_k]$ for $1 \le k \le 10$. > > (b) Find $\text{Cov} (X_k,X_5)$ for $1 \le k \le 10$. Trial:I think $X_k\sim \text{Bin}(n=100,p=1/5)$ and $E(X_k)=np=20,V(X_k)=npq=16$.Am I right? Lastly how I find $\text{Cov}(X_k,X_5)$.


title: Unit root test specifications?

text: I am puzzled as to what specification I should include in my unit root tests of the following data: ![enter image description here](http://i.stack.imgur.com/kXSXg.jpg). I will use ADF, KPSS and ZA tests. I can see there is a break in trend at observation 9. However, for the former two tests, especially for ADF, I do not know what to choose as a final specification.


title: Non-nested model verification

text: My questions are two fold:   1. Is there a generally accepted statistic used to compare non-nested, nonlinear models with different numbers of parameters? I'm thinking RMSE, but wondering what other options might be?   2. I have a small dataset and am hoping to classify curves from it into three different families (whose curves I know ahead of time). I am hoping to use Leave One Out cross-validation to classify each curve, but am worried that the variance of GOF statistics will be too large. What are some other options? Thanks so much!


title: How to derive a confidence interval from an F distribution?

text: So, this is the question I'm working on: Suppose we observe a random sample of five measurements: 10, 13, 15, 15, 17, from a normal distribution with unknown mean $µ_1$ and unknown variance $σ_1^2$. A second random sample from another normal population with unknown mean $µ_2$ and unknown variance $σ_2^2$ yields the measurements: 13, 7, 9, 15, 11. b. Use the pivotal method (and a pivotal statistic with F distribution) to derive a 95% confidence interval for $σ_2/σ_1$. Work it out for these data. And test the null hypothesis that $σ_2 = σ_1$ at the 5% level of significance. [6] (recall that $F_{\nu_1,\nu_2,a} = 1/F_{\nu_2,\nu_1,1-a}$). So, I'm completely at a loss as to how I can use the pivotal method on the F distribution. Please help me.


title: Blocked repeated measures ANOVA in R

text: I am trying to test if richness varies by treatment (Severity) over time (Block) using a repeated measures ANOVA in R. Any suggestions on how to do this correctly? I have tried, but get an error message stating I have an extra = in the formula:               TWP.aov <- aov(Richness ~ Severity * Year + Error(Plot/Severity), data = TWP)     summary(TWP.aov)       Here is a subset of the data (called TWP which I read in as a csv file):               Block       Severity    Plot    Richness     2003-2004   High        A       18     2003-2004   High        B       24     2003-2004   High        C       21     2005-2006   High        A       28     2005-2006   High        B       24     2005-2006   High        C       20     2007-2009   High        A       14     2007-2009   High        B       27     2007-2009   High        C       29     2003-2004   Low         A       12     2003-2004   Low         B       10     2003-2004   Low         C       14     2005-2006   Low         A       18     2005-2006   Low         B       16     2005-2006   Low         C       14     2007-2009   Low         A       8     2007-2009   Low         B       19     2007-2009   Low         C       20      When I input the above subset I get the correct summary table:               Error: Plot               Df Sum Sq Mean Sq F value Pr(>F)     Residuals  2  49.33   24.67                  Error: Plot:Severity               Df Sum Sq Mean Sq F value Pr(>F)       Severity   1 304.22  304.22   85.56 0.0115 *     Residuals  2   7.11    3.56                      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Error: Within                    Df Sum Sq Mean Sq F value Pr(>F)     Block           2  43.00  21.500   0.745  0.505     Severity:Block  2   1.44   0.722   0.025  0.975     Residuals       8 230.89  28.861        My actual data's summary table looks like this though:               Error: Plot                    Df Sum Sq Mean Sq F value Pr(>F)       Severity        1    248  247.66   4.013 0.0479 *     Block           2    493  246.67   3.997 0.0214 *     Severity:Block  2    244  122.17   1.980 0.1436       Residuals      99   6110   61.71                      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      As of June 12, 2013: As an update, I still have not resolved the "Error() model is singular" issue with the aov function, so I have tried to run this in the nlme package using the following code:               TWP.lme <- lme(Richness ~Severity * Block, random = ~1|(Plot/Severity), data = TWP)         summary(TWP.lme)         anova(TWP.lme)      I do not get any error message, and I believe I have structured the error term correctly. I have checked the results visually against the graphs of the data and all appears correct. Is there a way to check the results are correct? Why is there no error message regarding the model is singular?


title: How to do a Bayesian survival analysis and determine which variables are useful?

text: I can measure two variables on each patient in a survival study (I have the measurements and the survival times; some patients outlive the study and are therefore censored). I know that it is possible to use either of these variables to assign patients to groups with different survival curves (using a simple threshold on the variables). What I want to know is whether using both variables allows me to assign patients better than just one of the variables — i.e., whether the survival characteristics are more different when both variables are included; whether a given variable provides information that the other one doesn't. I would prefer to do this in a Bayesian manner using JAGS (BUGS); I tend to find these models easier to understand. Thanks in advance.


title: Activation value at output neuron equals 1, and the network doesn't learn anything

text: I'm implementing a typical neural network with 1 hidden layer. The network does well with the logic XOR and other simple problems, but fails miserably when encountering a (16-input, 20~30 hidden, 3 output) training set. The problem is that the network does not learn anything, even with learning rate = 1.0 and momentum 0.8 ~ 0.9. I did a lot of debug and found that at some output, the activation value is 1.0, even though mathematically the result of the logistic function (1/(1+exp(-signal)) can not be 1, but due to rounding, the computer produces 1.0 for signal that is >40. Now, with the activation value equals 1.0, the gradient is always 0 regardless of the desired output because the term `(1-activation)` will be 0, so in the equation:               delta = (1-activation) * (activation) * (activation - desired_output)      delta will be zero all the time, so no learning at all. I tried to initialize the weights at very small value (0.00001 or so) but eventually all of them reach 3.0 ~4.0, or even some extreme value like 30~40, so with 20~30 hidden neurons it's very easy to end up with output signal >40. Obviously I must have missed something, but with this kind of neural net (with hundreds of weights) I can't debug it too see what's wrong, so my question is: do all of the above symptoms signal anything wrong with my code? Here are some important code:               class NNWeight {         double ** weights;         double *bias;     }     double sigmoid(double value) {         return 1.0/(1.0+exp(-value));         }     double Dsigmoid(double value){  //gradient         return value*(1-value);     }         //wInpHid = input-to-hidden layer weights (of type NNWeight)         //wHidOut = hidden-to-output layer weights (of type NNWeight)              void updateWeight() {         for (int iInp = 0; iInp<nInp; iInp++)             for (int iHid = 0; iHid<nHid; iHid++)                 wInpHid.weights[iInp][iHid] += dInpHid.weights[iInp][iHid];         for (int iHid = 0; iHid<nHid; iHid++) {             wInpHid.bias[iHid] += dInpHid.bias[iHid];         }         for (int iHid = 0; iHid<nHid; iHid++)             for (int iOut = 0; iOut<nOut; iOut++)                 wHidOut.weights[iHid][iOut] += dHidOut.weights[iHid][iOut];         for (int iOut = 0; iOut<nOut; iOut++)             wHidOut.bias[iOut] += dHidOut.bias[iOut];     }     void feedFW() {         bool correct = true;         // Feed the hidden layer         for (int iHid = 0; iHid<nHid; iHid++) {             hid[iHid] = wInpHid.bias[iHid];             for (int iInp = 0; iInp<nInp; iInp++) {                 hid[iHid] += inp[iInp] * wInpHid.weights[iInp][iHid];             }             hid[iHid] = sigmoid(hid[iHid]);         }         //Feed the output layer         for (int iOut = 0; iOut<nOut; iOut++){             out[iOut] = wHidOut.bias[iOut];             for(int iHid = 0; iHid < nHid; iHid++)                 out[iOut] += hid[iHid]*wHidOut.weights[iHid][iOut];                  out[iOut] = sigmoid(out[iOut]);         }     }     void backProp() {         double SE =0.0;         //Calc the output errors:         for (int iOut = 0; iOut<nOut; iOut++){             SE += (out[iOut] - tar[iOut])*(out[iOut] - tar[iOut]);             outErr[iOut] = (tar[iOut] - out[iOut]) * Dsigmoid(out[iOut]);         }         //Calc the hidden errors:         for (int iHid = 0; iHid<nHid; iHid++){             double sum_wHidOut = 0.0;             for (int iOut = 0; iOut<nOut; iOut++) {                 sum_wHidOut += outErr[iOut]*wHidOut.weights[iHid][iOut];             }             hidErr[iHid] = sum_wHidOut * Dsigmoid(hid[iHid]);         }         //Calc hidden - output layer delta (dHidOut)         for (int iHid = 0; iHid<nHid; iHid++)             for(int iOut = 0; iOut<nOut; iOut++){                 dHidOut.weights[iHid][iOut] = lRate * outErr[iOut] * hid[iHid] + dHidOut.weights[iHid][iOut] * momentum; //lRate = learning rate             }         for (int iOut = 0; iOut<nOut; iOut++){             dHidOut.bias[iOut] = lRate * outErr[iOut] + dHidOut.bias[iOut] * momentum;         }         //Calc input - hidden layer delta (dInpHid)         for (int iInp = 0; iInp<nInp; iInp++)             for (int iHid = 0; iHid<nHid; iHid++)                 dInpHid.weights[iInp][iHid] = lRate * hidErr[iHid] * inp[iInp] + dInpHid.weights[iInp][iHid] * momentum;         for (int iHid = 0; iHid<nHid; iHid++)             dInpHid.bias[iHid] = lRate * hidErr[iHid] + dInpHid.bias[iHid] * momentum;                 //Update the weights         updateWeight();     }


title: Zero values and discontinuity in explanatory variable

text: One of my independent variables measures worker productivity through the variable $\frac{\log{sales}}{\text{# of workers}}$, and I'm creating one variable for skilled and another for unskilled workers. Past a point, as the number of skilled workers decreases there is a strong correlation with an increase in productivity per skilled worker, which is why I'm thinking it would create a bias if I simply let productivity be zero when the number of workers is zero. At one worker the variable has a very high value and at zero workers, it would be equal to zero. At the same time, throwing observations with zero workers of that category out is not an acceptable solution. I'm using a GLM with an identity link function. How can I account for these observations (perhaps in a separate dummy variable?) without biasing the coefficient?


title: Difference between a stochastic and deterministic optimization problem?

text: I am reading a lot about classifying optimization problems. But i cannot find a real life example that treads the difference between a stochastic and deterministic optimization problem.


title: CDF of conditional distribution

text: Suppose that $X$ and $Y$ are iid normally distributed and $a$ is a scalar. What is $\Pr(Y+aX<0 | X>0)$?


title: Testing for moderation with continuous vs. categorical moderators

text: I am testing an interaction effect where $X$ and $Y$ are continuous variable and $M$ (Moderator) is a categorical variable (effects coding $+1$, $-1$). I have no clue about how to do a post-hoc probing of slopes in this case. For a continuous moderator we calculate $Z_{\rm above}$ and $Z_{\rm low}$, calculate their crossproduct with $X$ (e.g., $Z_{\rm Above}*X$ and $Z_{\rm low}*X$) and then run the regression with their combinations of ${\rm High}$ and ${\rm Low}$ with $Y$. I do not see that the same method above could be applied for a categorical moderator. How can I test that?


title: Negative values for AIC in General Mixed Model

text: I'm trying to select the best model by the AIC in the General Mixed Model test. The best model is the model with the lowest AIC, but all my AIC's are negative!   * So is the biggest negative AIC the lowest value?   * Or is the smallest negative AIC the lowest value, because it's closer to 0? For example is AIC -201,928 or AIC -237,847 the lowest value and thus the best model?


title: What is the correct test to check whether a observable has changed after modifying a condition?

text: I run the following 'experiment':   * take measurement x of observable   * change one parameter   * take measurement y of observable again I do this many times and want to test whether changing the parameter has an effect on my observables. More specifically, I want to test whether changing the parameter does increase the measured value, that is, is `x < y`? As far as I understand correctly, what I should do is a paired t-test, since I'm comparing the same measurement, before and after I change something. However, what I'm unsure about, is what statistical test is the correct one to use?               [h, p_ttest1] = ttest(all_before, all_after, [] , 'left');     [h, p_ttest2] = ttest2(all_before, all_after, [], 'left');     [h, p_ttest3] = ttest(all_after - all_before, [] , 'right');     [p_signrank, h] = signrank(x, y);       the left-sided ttest and ttest2 yield very different p values, that's why I'm trying to understand what test is the correct one to use here.


title: The exogeneity condition and correlation

text: My book says that the exogneity condition, the condition that should hold in order for $b$ to be consistent, is the following: $$\operatorname{plim}(\frac{1}{n}X'\epsilon)=0.$$ It then says that the $j$th component of this vector can be written as $\frac{1}{n} \sum_{i=1}^n x_{ji}\epsilon_i$, so that the condition is equivalent to that the explanatory variables should be asymptotically uncorrelated with the disturbances. I don't however see why $\frac{1}{n} \sum_{i=1}^n x_{ji}\epsilon_i$ represents the correlation between the $x$ and \epsilon$ variables. Could anyone please explain me why this is so?


title: Gaussian Ratio Distribution: Derivatives wrt underlying $\mu$'s and $\sigma^2$s

text: I'm working with two independent normal distributions $X$ and $Y$, with means $\mu_x$ and $\mu_y$ and variances $\sigma^2_x$ and $\sigma^2_y$. I'm interested in the distribution of their ratio $Z=X/Y$. Neither $X$ nor $Y$ has a mean of zero, so $Z$ is not distributed as a Cauchy. I need to find the CDF of $Z$, and then take the derivative of the CDF with respect to $\mu_x$, $\mu_y$, $\sigma^2_x$ and $\sigma^2_y$. Does anyone know a paper where these have already been calculated? Or how to do this myself? I found the formula for the CDF in a 1969 paper, but taking these derivatives will definitely be a huge pain. Maybe someone has already done it or knows how to do it easily? I mainly need to know the signs of these derivatives. This paper also contains an analytically simpler approximation if $Y$ is mostly positive. I can't have that restriction. However, maybe the approximation has the same sign as the true derivative even outside the parameter range?


title: Incorporating a response variable in principal component analysis

text: I'm a biologist, and I have a large dataset that I'm trying to analyze. Here are the variables I'm working with:   * levels of 211 different metabolites in 16 different blood samples (predictor variables)   * how well each of the 16 blood samples performed in a specific test (response variable) I am trying to figure out which variable(s) are the most important in predicting the performance of a blood sample in the test. I would like to make these data more manageable by doing a PCA, but I'm new to this sort of analysis. I understand that a PCA will create groups of principal components (it will group metabolites that covary with each other and label each of these groups a principal component), but I'm not sure how to take into account the response variable in this analysis. Any help would be much appreciated!!! Thank you.               > summary(metabolites_princomp)     Importance of components:                               Comp.1     Comp.2      Comp.3      Comp.4     Standard deviation     3.9608225 0.40128486 0.259868774 0.215004349     Proportion of Variance 0.9805072 0.01006435 0.004220736 0.002889179     Cumulative Proportion  0.9805072 0.99057153 0.994792267 0.997681446                                Comp.5       Comp.6       Comp.7     Standard deviation     0.12944496 0.0835343589 0.0745499634     Proportion of Variance 0.00104725 0.0004361243 0.0003473561     Cumulative Proportion  0.99872870 0.9991648202 0.9995121763                                  Comp.8       Comp.9      Comp.10     Standard deviation     0.0606023885 0.0467885157 2.566109e-02     Proportion of Variance 0.0002295406 0.0001368228 4.115572e-05     Cumulative Proportion  0.9997417169 0.9998785397 9.999197e-01                                 Comp.11      Comp.12      Comp.13     Standard deviation     2.203313e-02 1.872209e-02 1.450158e-02     Proportion of Variance 3.034116e-05 2.190729e-05 1.314349e-05     Cumulative Proportion  9.999500e-01 9.999719e-01 9.999851e-01                                 Comp.14      Comp.15      Comp.16     Standard deviation     1.141596e-02 7.517519e-03 7.194785e-03     Proportion of Variance 8.145258e-06 3.532068e-06 3.235308e-06     Cumulative Proportion  9.999932e-01 9.999968e-01 1.000000e+00


title: Repeated measures ANOVA, unbalanced and missing values in R

text: I have chemical compositions of individuals collected in 6 different sites in monthly intervals over a period of 6 months. I would like to determine whether there is significant variation among sites as well as among months within sites. I was advised to use a repeated measures ANOVA to address this. However my design is unbalanced and has missing values. So, for example not all the sites are represented at all times and the number of individuals per site and time also varies. Is it possible to do a repeated measures ANOVA with this dataset? Should I use the lmer function from the lme4 package?


title: Alternatives to Hypothesis Testing in R

text: I've been reading a number of articles and blog posts about issues relating to hypothesis testing. While these sources seem to put forth legit problems with hypothesis testing and it's interpretations, I'm not seeing anything in terms of alternatives to it (beyond bayesian hypothesis testing). What are some alternative to hypothesis testing in R? Problems with the hypothesis testing approach So-called Bayesian hypothesis testing is just as bad as regular hypothesis testing Alternatives to Statistical Hypothesis Testing EDIT: For example, let's consider the following situation. We have monthly data for one year where conversions_a is the control and conversions_b is the experimental data for conversion when using a different headline.               df = data.frame(year_one_a = 1:12, conversions_a = rnorm(12),                      year_one_b = 1:12, conversions_b = rnorm(12)+5)


title: Interpreting odds ratios

text: I have the following set of results for one of the factors (birth weight) with different levels and their corresponding Odds ratios for survival. I am using the first level (<1.25) as the reference level: Birth weight (kg):               Levels      Number/level    Odds Ratio     <1.25       1615            1.00     1.25-1.46   1617            1.37     1.46-1.65   1462            1.25     1.65-1.87   1632            1.68     >1.87       1466            2.35      From this result, I am trying to estimate, by using the OR, the approximate number of newborns that would survive in each level. Are there ways that this can be achieved using R? In a separate calculation, the actual number of survival per level are shown below:               Levels      Number/level    Odds Ratio    Actual number survived     <1.25       1615            1.00         1088     1.25-1.46   1617            1.37         1346     1.46-1.65   1462            1.25         1238     1.65-1.87   1632            1.68         1447     >1.87       1466            2.35         1351      EDIT: The above is just one of the independent factors from my model.               mod <- glm(formula = surv ~ as.factor(var1) + as.factor(var2)+...as.factor(varn)+family = binomial(link = "logit"), data=mydf)               summary(mod)          glm(formula = surv ~ as.factor(season) + as.factor(bwt5) + as.factor(prectem5) +          as.factor(pcscore) + as.factor(pindx5) + as.factor(presp2) +          as.factor(ppscore) + as.factor(mtone2) + as.factor(fos) +          as.factor(psex) + as.factor(pscolor) + as.factor(pshiv) +          as.factor(backfat5) + as.factor(srect2) + as.factor(gest3) +          as.factor(int3) + as.factor(agit) + as.factor(tacc), family = binomial(link = "logit"),          data = lesna)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -3.0562   0.3120   0.4478   0.5929   1.9412            Coefficients:                           Estimate Std. Error z value Pr(>|z|)         (Intercept)           1.441842   0.290536   4.963 6.95e-07 ***     as.factor(season)2   -1.064053   0.107666  -9.883  < 2e-16 ***     as.factor(bwt5)2      0.314332   0.099776   3.150 0.001631 **      as.factor(bwt5)3      0.223824   0.110566   2.024 0.042935 *       as.factor(bwt5)4      0.524586   0.120182   4.365 1.27e-05 ***     as.factor(bwt5)5      0.854196   0.138993   6.146 7.97e-10 ***     as.factor(prectem5)2  0.745025   0.094238   7.906 2.66e-15 ***     as.factor(prectem5)3  0.856777   0.098326   8.714  < 2e-16 ***     as.factor(prectem5)4  0.997219   0.111529   8.941  < 2e-16 ***     as.factor(prectem5)5  0.930925   0.120052   7.754 8.88e-15 ***     as.factor(pcscore)2   0.384534   0.137564   2.795 0.005185 **      as.factor(pcscore)3   0.668390   0.154608   4.323 1.54e-05 ***     as.factor(pindx5)2    0.243485   0.101755   2.393 0.016718 *       as.factor(pindx5)3    0.262779   0.108809   2.415 0.015733 *       as.factor(pindx5)4    0.595672   0.118124   5.043 4.59e-07 ***     as.factor(pindx5)5    0.467277   0.120401   3.881 0.000104 ***     as.factor(presp2)2   -0.286214   0.126012  -2.271 0.023127 *       as.factor(ppscore)2  -0.246369   0.093568  -2.633 0.008462 **      as.factor(mtone2)2   -0.482397   0.118218  -4.081 4.49e-05 ***     as.factor(fos)2      -0.255652   0.075749  -3.375 0.000738 ***     as.factor(psex)2      0.182437   0.066964   2.724 0.006442 **      as.factor(pscolor)2  -0.694197   0.282069  -2.461 0.013852 *       as.factor(pshiv)2    -0.241515   0.080792  -2.989 0.002796 **      as.factor(backfat5)2 -0.309427   0.104176  -2.970 0.002976 **      as.factor(backfat5)3 -0.004152   0.108669  -0.038 0.969523         as.factor(backfat5)4  0.013233   0.103491   0.128 0.898257         as.factor(backfat5)5 -0.221935   0.104639  -2.121 0.033926 *       as.factor(srect2)2   -0.236981   0.104962  -2.258 0.023960 *       as.factor(gest3)2     0.207375   0.106065   1.955 0.050562 .       as.factor(gest3)3     0.904959   0.191307   4.730 2.24e-06 ***     as.factor(int3)2     -0.204870   0.127674  -1.605 0.108573         as.factor(int3)3     -1.271092   0.388924  -3.268 0.001082 **      as.factor(agit)2     -0.496856   0.157553  -3.154 0.001613 **      as.factor(agit)3     -0.360247   0.148520  -2.426 0.015284 *       as.factor(tacc)2     -0.282180   0.090556  -3.116 0.001833 **      as.factor(tacc)3     -0.429249   0.082520  -5.202 1.97e-07 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for binomial family taken to be 1)              Null deviance: 7096.2  on 7791  degrees of freedom     Residual deviance: 6143.0  on 7756  degrees of freedom     AIC: 6215          Number of Fisher Scoring iterations: 5          exp(mod$coefficients) # odds ratios              (Intercept)   as.factor(season)2     as.factor(bwt5)2                 4.2284770            0.3450544            1.3693441          as.factor(bwt5)3     as.factor(bwt5)4     as.factor(bwt5)5                 1.2508506            1.6897596            2.3494840      as.factor(prectem5)2 as.factor(prectem5)3 as.factor(prectem5)4                 2.1064944            2.3555559            2.7107340      as.factor(prectem5)5  as.factor(pcscore)2  as.factor(pcscore)3                 2.5368547            1.4689297            1.9510940        as.factor(pindx5)2   as.factor(pindx5)3   as.factor(pindx5)4                 1.2756866            1.3005395            1.8142499        as.factor(pindx5)5   as.factor(presp2)2  as.factor(ppscore)2                 1.5956440            0.7511016            0.7816335        as.factor(mtone2)2      as.factor(fos)2     as.factor(psex)2                 0.6173022            0.7744113            1.2001382       as.factor(pscolor)2    as.factor(pshiv)2 as.factor(backfat5)2                 0.4994756            0.7854368            0.7338673      as.factor(backfat5)3 as.factor(backfat5)4 as.factor(backfat5)5                 0.9958568            1.0133207            0.8009678        as.factor(srect2)2    as.factor(gest3)2    as.factor(gest3)3                 0.7890065            1.2304445            2.4718305          as.factor(int3)2     as.factor(int3)3     as.factor(agit)2                 0.8147530            0.2805250            0.6084405          as.factor(agit)3     as.factor(tacc)2     as.factor(tacc)3                 0.6975039            0.7541382            0.6509975           exp(coef(mod)) #exponentiated coefficients     exp(confint(mod)) # 95% CI      RUNNING THE SCRIPT WITHOUT type="response"               > (pred1s <- predict(mod,newdata=as.data.frame(with(lesna,list(season=1,bwt5=1,     +            prectem5=1,pcscore=1,pindx5=1,presp2=1,ppscore=1,mtone2=1 .... [TRUNCATED]             1      1.441842      RUNNING THE SCRIPT WITH type="response":               > (pred1s <- predict(mod,newdata=as.data.frame(with(lesna,list(season=1,bwt5=1,     +            prectem5=1,pcscore=1,pindx5=1,presp2=1,ppscore=1,mtone2=1 .... [TRUNCATED]              1      0.8087397


title: Methods for sampling mutual independent and k-wise independent random variables

text: 1. I wonder how a sequence of mutually independent random variables (each with a different or same distribution) is sampled? In other words, how can a sequence of values be ensured to be samples of a sequence of independent random variables?   2. Furthermore how is a sequence of k-wise independent random variables is sampled? I would like to know the basic methods, and what sources for that purpose and going beyond that. Thanks!


title: How to conduct predictor selection in a generalized linear mixed model?

text: I have 18 predictors in a binary generalized linear mixed model (repeated measurements, over a 1000 subjects). I would like to trim the model a bit and remove some noise and useless predictors. Unfortunately, `PROC GLIMMIX` does not have any facility to do this. I could not find an R package that would do this (`step()` function style). If I were to try this manually, say begin with a full model (all predictors in) and do a 'backward selection', what criterion could I use to do this quickly? Could I use, say, p-values? But at what significance level?


title: Why am I getting a zero parameter estimate and SE for one of my variables in a logistic regression?

text: I am doing a binary logit analysis, in which I'm trying to fit a model to my data to explain why some towns adopt open space subdivision ordinances (OP), using a handful of discrete and continuous independent variables. One of the explanatory variables, median household income, is coming back with an estimated beta of 0, with a S.E. of 0, and a p value of (0.14), and the Odds Ratio is 1 (S.E. of 0). It is highly unlikely that there is absolutely no relation between household income and towns passing a growth management ordinance, so this seems pretty odd. The rest of my variables have normal looking estimates and errors. This is my first time asking a question on this site, so please forgive me if I leave out critical information - I'll explain more if needed. I tried to cut and paste my output below, but the formatting of the table disappeared, making it impossible to decipher. Also, my descriptive statistics for household income look fine, with realistic minimum, maximum and mean numbers, and the original data doesn't show anything obviously wrong with the data itself, which makes me think I'm doing something wrong when estimating my model. On a side note, my odds ratio results for percent change in housing from 1990-2000 (PCHU90_00) look really high, and well, odd =) I'm not sure what's going on there, either. Any suggestions would be appreciated. Thank you! Here's the model:               LOGIT     MODEL OP = CONSTANT+COAST+LA_SQMI+PCHU00_10+PCHU70_80+PCHU80_90+PCHU90_00+CONSTEMP+          HHINCOME+LT+PCHPOP00_10+PCHPOP60_70+PCHPOP70_80+PCHPOP80_90+PCHPOP90_00+PCTDEM      And the results for the two variables I'm concerned about:               HHINCOME: 0 Beta (0 SE), z=1.46, p=0.14, 1.0 Odds Ratio (0 SE)     PCHU90_00: 9.64 (3.89 SE), z=2.47, p=.01, 15310.11 Odds Ratio (59625.59 SE)      Model fit info seems okay -               Log-Likelihood of Constants only Model = LL(0):-83.14     2*[LL(N)-LL(0)] :37.13     df:15     p-value:0.00     McFadden's Rho-squared  0.22     Cox and Snell R-square  0.25     Naglekerke's R-square   0.35


title: R packages or open source software for training Hidden Markov chains

text: Are there any well-designed R packages or other open-source software for training Hidden Markov chains?


title: What can I conclude from the following one way ANOVA result?

text: I used a one-way ANOVA to test if importance levels of attributes have changed over the years (Data comes from a survey of consumers) with the following results (Mean ratings):               +-----+---------+     | Year|Attribute|     +-----+---------+     | 2000|     4.39|     +-----+---------+     | 2010|     4.30|     +-----+---------+     | 2011|     4.38|     +-----+---------+     | 2012|     4.43|     +-----+---------+      The overall ANOVA is significant and the Tukey post-hoc test of pairwise comparisons indicates that the mean rating for the year 2012 is greater than the mean rating for the year 2010 at $\alpha=0.05$ (i.e., $4.43 > 4.30$). None of the other comparison pairs are significant. How can I translate the ANOVA results into plain English? It seems an odd result. Importance levels seem stable for the years 2009, 2010 and 2011. But, the importance level for the year is greater for the year 2012 than the importance level for the year 2010 but not for the year 2009 and 2011. I understand the statistics as to why this can happen but am confused with how to translate the above into plain English.


title: Estimating credit card default probabilities

text: Folks, I am working on a credit card defaults and transition probabilities. For example a single credit card account could be in a number of states: up- to-date, 30, 60, 90 days in arrears or in default.   * Are there packages in R that do estimation of the transition probabilities given historical monthly cohort defaults?   * Any pointers to papers specific to this type of estimation?   * Simulation of future "paths" of defaults. Note that transitions occur discretely on a monthly basis. Thanks much for your time, KW