title: Finding marginal densities of $f (x,y) = c \sqrt{1 - x^2 - y^2}, x^2 + y^2 \leq 1$

text: As the title says, I'm looking for the marginal densities of $$f (x,y) = c \sqrt{1 - x^2 - y^2}, x^2 + y^2 \leq 1.$$ So far I have found $c$ to be $\frac{3}{2 \pi}$. I figured that out through converting $f(x,y)$ into polar coordinates and integrating over $drd\theta$, which is why I'm stuck on the marginal densities portion. I know that $f_x(x) = \int_{-\infty}^\infty f(x,y)dy$, but I'm not sure how to solve that without getting a big messy integral, and I know the answer isn't supposed to be a big messy integral. Is it possible to instead find $F(x,y)$, and then take $\frac{dF}{dx}$ to find $f_x(x)$? That seems like the intuitive way to do it but I can't seem to find anything in my textbook that states those relationships, so I didn't want to make the wrong assumptions.


title: How to name the ticks in a python matplotlib boxplot

text: Python matplotlib has a boxplot command. Normally, all the parts of the graph are numerically ticked. How can I change the ticks to names instead of positions? For illustration, I mean the Mon Tue Wed labels like in this boxplot: ![alt text](http://i.stack.imgur.com/ZFENf.png)


title: Method for estimating (thermodynamic) scaling fit parameters

text: Let's say I have some measurements, label them $N$, which depend on the variables $x$ and $y$. Further, I hypothesize that there is a scaling relationship between $N$, $x$, and $y$ of the form $x^{\alpha}N = f(x^{\beta}(y-\gamma))$, where $\alpha$, $\beta$, and $\gamma$ are unknown parameters and $f$ is some unknown function. If I have many values of $N$, recorded for many values of $x$ and $y$, is there a way to estimate the values of the parameters $\alpha$, $\beta$, and $\gamma$? I don't really care what the function $f$ is. I'm sorry if this is unclear, or simplistic; I'm embarrassed to admit I know very little statistics. It would be extremely helpful if you know of an existing computational library (e.g. an R package) that would be helpful here, too. Thanks!


title: Load Linux workspace in S-PLUS for Windows

text: Is it possible to load an S-PLUS Linux workspace in Windows? If I try it I get this error: "Problem in exists(name, where = db): This directory has both Unix style __ _nonfile and Windows style __ _nonfi" The __nonfi file is created when I first try to load that Linux workspace in Windows. Is there any way to convert it to a Windows workspace? I'm using S-PLUS 8.0.


title: phantom NA/NaN/Inf in foreign function call (or something altogether different?)

text: Please forgive the puzzled title and the length of this message - I thought it would be best to be as complete as possible and to show the avenues I have explored. I'm trying to fit a linear model to data with a binary dependent variable (i.e. Target.ACC: accuracy of response) using lrm, and thought I would start from the most complex model (of which "sample1.lrm1" is a trimmed version). I got the error shown below. (sample1 is available at http://tinyurl.com/bwqq7ya) For info:               str(sample1)     >'data.frame':   14022 obs. of  5 variables:     > $ Target.ACC : Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...     > $ Word.Order : Factor w/ 2 levels "HeadMod*","ModHead": 1 1 1 1 1 1 1 1 1 1 ...     > $ Target.RESP: Factor w/ 2 levels "1","2": 2 2 2 2 2 2 2 2 2 2 ...     > $ L1         : Factor w/ 3 levels "English","German",..: 2 2 2 2 2 2 2 2 2 2 ...     > $ Relation   : Factor w/ 4 levels "For","From","MadeOf",..: 1 1 1 1 1 1 1 1 1 1 ...      Commands and error message:               library(rms)     sample1.dd = datadist(sample1)     options(datadist="sample1.dd")     sample1.lrm = lrm(Target.ACC ~ (L1 + Relation + Target.RESP + Word.Order)^2,                        sample1, x=T, y=T)     > Error in lrm(Target.ACC ~ (L1 + Relation + Target.RESP + Word.Order)^2,  :     > Unable to fit model using “lrm.fit”      So I tried to narrow down the error by looking at all the combinations manually, and the problem appears to be specifically with the interaction between Word.Order and Target.RESP. Models including interaction of these variables with other variables (e.g. L1, Relation) can be fitted without problem.               sample1.lrm = lrm(Target.ACC ~ (Target.RESP + Word.Order)^2, sample1, x=T, y=T)     > Error in lrm(Target.ACC ~ (Target.RESP + Word.Order)^2, dat, x = T, y = T) :     > Unable to fit model using “lrm.fit”      unproblematic:                 sample1.lrm1 = lrm(Target.ACC ~ (L1 + Relation + Target.RESP)^2, sample1, x=T,y=T)       sample1.lrm2 = lrm(Target.ACC ~ (L1 + Relation + Word.Order)^2, sample1, x=T,y=T)      When running the problematic analysis on a smaller sample of the same data, I get a different (more precise?) error message:               sample2 <- sample1[1:500,]     sample2.lrm = lrm(Target.ACC ~ (Target.RESP + Word.Order)^2, sample2, x=T, y=T)     > Error in fitter(X, Y, penalty.matrix = penalty.matrix, tol = tol, weights = weights,:     > NA/NaN/Inf in foreign function call (arg 1)      But I cannot find any NA in the data:               table(complete.cases(sample2))     > TRUE     > 500      Some portions of the data don't appear to contain any of the offending "bit":               sample3 <- sample1[12500:13000,]     sample3.lrm = lrm(Target.ACC ~ (Target.RESP + Word.Order)^2, sample3, x=T, y=T)      Could one of you shine your light on this puzzle, please? If that includes pointing me towards some background reading, that would be great too. Many thanks in advance.               sessionInfo()     > R version 2.15.0 (2012-03-30)     > Platform: i386-pc-mingw32/i386 (32-bit)     > locale:     > [1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United Kingdom.1252        > [3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C                                > [5] LC_TIME=English_United Kingdom.1252         > attached base packages:     > [1] splines   stats     graphics  grDevices utils     datasets  methods   base          > other attached packages:     > [1] rms_3.5-0        Hmisc_3.9-3      survival_2.36-12     > loaded via a namespace (and not attached):     > [1] cluster_1.14.2 grid_2.15.0    lattice_0.20-6 tools_2.15.0


title: Random generation of ARMA(2,2) Gaussian time series

text: I get very poor replication of longitudinal parameters from my own program using the Box-Jenkins model. I had no such problem with my own program generating AR(1) Gaussian data. Is there some trick that I haven't discovered yet? I would like to use the validation framework of Martinez-Rivera & Ventosa-Santaularia (2012) to check a new time series regression method that is better than the usual on AR(1), but itself depends on ARMA random generation when applied to ARMA data.


title: Drawing from a conditional density

text: I have a simple question. Suppose $X=(X_1,X_2,X_3)$ is multivariate normal. What's the best (quickest) way to draw from the conditional density $X_1\mid \exp(X_1)+\exp(X_2)+X_3$?


title: Can you recommend a book to read before Elements of Statistical Learning?

text: Based on this post, http://quant.stackexchange.com/questions/111/how-can-i-go- about-applying-machine-learning-algorithms-to-stock-markets, I want to digest Elements of Statistical Learning. Fortunately it is available for free and I started reading it. I don't have enough knowledge to understand it. Can you recommend a book that is a better introduction to the topics in the book? Hopefully something that will give me the knowledge needed to understand it? ### Related > How to get an introductory understanding of Machine Learning?


title: Are support vector regression and kernel ridge regression used for the same type of problems?

text: Does, in some sense, Kernel Ridge Regression refer to the same class of problems as Support Vector Regression? Should I use them almost equivalently?


title: Cluster Analysis followed by Discriminant Analysis

text: What is the rationale, if any, to use Discriminant Analysis (DA) on the results of a clustering algorithm like k-means, as I see it from time to time in the literature (essentially on clinical subtyping of mental disorders)? It is generally not recommended to test for group differences on the variables that were used during cluster construction since they support the maximisation (resp. minimisation) of between-class (resp. within-class) inertia. So, I am not sure to fully appreciate the added value of predictive DA, unless we seek to embed individuals in a factorial space of lower dimension and get an idea of the "generalizability" of such a partition. But even in this case, cluster analysis remains fundamentally an exploratory tool, so using class membership computed this way to further derive a scoring rule seems strange at first sight. Any recommendations, ideas or pointers to relevant papers?


title: Forecasting a "chaotic" time series

text: Here are four graphs, 1, autocorrelation, autocovariance, partial-correlation and cross-correlation calculated from a time series are given. 2, The time series I need to do some predictions on them. My questions are 4 : 1, I do not quite understand what these four graphs can tell or show. Can anyone explain? 2, Does the four subplots tell if it is appropriate to do forecasting for the time series 3, If so, any method to do forecasting? 4, I am wondering if particle filter can be used. I am a little clueless about this problem. ![descriptive stats](http://i.stack.imgur.com/K18Aa.png) ![The time series which seems "chaotic"](http://i.stack.imgur.com/Gc1ZU.png)


title: Number of tests for rare events and multiple comparisons problem?

text: I am conducting some statistical tests on a rare event (11 cases out of 171 patient sample). I have two main hypotheses, but in the mean time I also want to conduct some exploratory tests on 6 other variables, in order to explore possible associations. (1) I am wondering if there are any rules about the number of tests that you can do based on the number of events. (2) Normally when we do one or several statistical tests, we consider the statistical power (which is related with sample size), and with several tests we also consider the multiple comparisons problem (which is related with the number of tests). In the case of rare events, the statistical power of each of my test would be low, but will rare event also influence (or exaggerate) the multiple comparisons issue, i.e. after conducted several tests, the chance of finding false positives are much higher for rare events?


title: Questions about acceptance rule of Metropolis algorithm

text: I have a question concerning about one step in the Metropolis algorithm. The algorithm proceeds as following,   1. Generate a proposed new sample value from the jumping distribution $Q(x'|x_t)$   2. Calculate the acceptance ratio $a = \frac{P(x')}{P(x_t)}$   3. If $a\geq 1$, accept $x'$ by setting $x_{t+1}= x'$   4. Else, accept $x'$ with probability $a$. That is, pick a uniformly distributed random number $r$ between 0 and 1, and if $r\leq a$ set $x_{t+1}=x'$, else set $x_{t+1} = x_t$ I feel unclear about step 4, in step 4 if we have $a<1$, which means $P(x_t)>P(x')$, $x'$ is further away from the large mass of the density than $x_t$. However, we are not dumping $x'$ directly, we leave it to a 'gamble' to decide its go or stay. This gamble is the uniform distribution, a random draw from uniform distribution then compare it with $a$. I do not understand what's happening here. Why won't we just dump $x'$, why we have to give it a chance to 'nature'?


title: Trend of a few time points

text: Suppose we are given a matrix of many rows (different genes for example) and few columns (different time points) and we want to identify the top rows (genes) that are following a trend, like monotonically increasing or decreasing. What is a good method to perform this in R? The answers providing also a statistical significance result (something like p-value) are more than welcome. Finally the methods to indicate more complicated cases - such as up-then-down treands - are really perfect.


title: How to calculate the weighted mean?

text: For example assume I have three observations as given in matrix A $A = \left[ \begin{array}{ccc} 1 & 0 & 1 \\ 1 & 1 & 0 \\ 0 & 1 & 1 \\ \end{array} \right]$ each row vector of matrix A has a weight value associated with it as $w = (0.5, 0.3, 0.2) $ Using the equation in Wikipedia which is: $$ \bar{x} = \sum_{i = 1}^{N} w_ix_i $$ I'm confused on the direction of $i$ in the above equation. I managed to interpret this in two ways below   1. $i$ represents a row in matrix A $$ \bar{x} = 0.5*\left[ \begin{array}{ccc} 1 & 0 & 1 \\ \end{array} \right] + 0.3*\left[ \begin{array}{ccc} 1 & 1 & 0\\ \end{array} \right] + 0.2*\left[ \begin{array}{ccc} 0 & 1 & 1 \\ \end{array} \right] $$ $$ \bar{x} = \left[ \begin{array}{ccc} 0.5 & 0 & 0.5 \\ \end{array} \right] + \left[ \begin{array}{ccc} 0.3 & 0.3 & 0\\ \end{array} \right] + \left[ \begin{array}{ccc} 0 & 0.2 & 0.2 \\ \end{array} \right] = \left[ \begin{array}{ccc} 0.8 & 0.5 & 0.7 \\ \end{array} \right] $$   2. $i$ represents a column in matrix A (column is represented as a row) $$ \bar{x} = \left[ \begin{array}{ccc} 0.5 & 0.3 & 0 \\ \end{array} \right] + \left[ \begin{array}{ccc} 0 & 0.3 & 0.2\\ \end{array} \right] + \left[ \begin{array}{ccc} 0.5 & 0 & 0.2 \\ \end{array} \right] $$ $$ \bar{x} = \left[ \begin{array}{ccc} 1.0 & 0.6 & 0.4 \\ \end{array} \right] $$ Which calculation is correct? For me the approach 1 seems correct but I'm not sure. Please explain the two scenarios   1. If the columns of matrix A represents attributes and the rows represent observations which method is correct?   2. If the columns of matrix A represents observations and the rows represent attributes which method is correct?


title: Confidence interval for the ratio of Poisson to normal variables

text: In an experiment, I'm counting failures on an semiconductor device for a given fluence (i.e. amount) of neutron radiation. The failures occur according to a Poisson law. The fluence is provided by the neutron facility. The fluence standard error $\Delta{x} = \frac{\text{s}}{\sqrt{n}}$ has been computed by the neutron facility during a beam calibration phase. $\text{s}$ is the sample standard deviation and $\text{n}$ (which I do not know) the number of calibration measurements done by the facility. We assume that the distribution of the measurements is normal. I would like to compute the 95% confidence interval for the failure rate $F$ which is the ratio of the number of failures over this fluence: $F = \text{number failures} / \text{fluence}$ When I disregarded the uncertainty on the fluence, I used the $\chi^2$ 95% bounds to provide lower and upper limits for the number of failures; but I'd like to integrate the uncertainty on the fluence as well. As an example, I could have 2 failures for a fluence of 1.2E6 neutrons. The fluence relative uncertainty ($\frac{\Delta{x}}{\bar{x}}$,computed from the standard error defined above) is given at 8% for this specific facility. What would be the 95% confidence interval for $F = 2 / 1.2\text{E}6$ ?


title: Maximizing returns - A Bayesian approach

text: I want to design a Bayesian model for a simple asset allocation problem. Say I can buy $a_i$ amounts of $N$ assets. The return values of these assets are given by random variables $r_i$ with known posteriors $p(r_i | \textbf{x})$, given some known data about the market $\textbf{x}$. I define my utility function (i.e. negative loss) as the total return of my portfolio: Total utility = $U(\textbf a, \textbf r) = \textbf{a} \cdot \textbf{r} = \sum_i a_i r_i$ where, sticking to the bayesian notation, $\textbf{a}$ are my actions, and $\textbf{r}$ are the random returns. I would like to maximize the total expected utility using bayesian decision theory. From what I understand, I could write this as follows: $\underset{\textbf{a}}{\operatorname{argmax}} E\{U(\textbf{a} ,\textbf{r}) | \textbf{x}\} = \underset{\textbf{a}}{\operatorname{argmax}} \int (\textbf{a} \cdot \textbf{r}) \; p(\textbf{r})\;d \textbf{r}$ Assuming that   1. $\sum_i a_i \le C$ is a constraint so that there is a limit in how many assets in total I can buy   2. I can only buy (not sell) assets, i.e. $a_i \ge 0$. My questions are:   1. How can I modify this formulation to **penalize risk** (uncertainty) in a principled manner ? (i.e. penalize the action of buying assets with high variance/entropy/uncertainty)   2. What type of optimization problem is this? What solvers are applicable to this problem? (assuming a finite number of assets $N$)   3. Are there readily available libraries for solving this problem? Perhaps in PyMC.? As a bonus extension:   1. How can I modify this formulation to add the fact that there is a limited random supply $q_i$ of asset $i$ that follows $p(q_i|\textbf{x})$?


title: Statistical thought experiment (possibly Bayesian) about survey sampling and propensity scores

text: For some practical application I recently came about the following thought experiment. Can anybody help? Suppose we administer a survey A to measure a variable $Y$. The response probabilities may depend on $Y$, as in $$P(R_a=1|Y)$$ where $R_a$ is the response indicator, 0 for nonresponse, 1 for response. We also administer a second survey B, for example using some other fieldwork protocol, which leads to $$P(R_b=1|Y).$$ Now assume we can estimate both response propensities, because we have $Y$ from some external source. The goal is now to find a sample which in expectation has response probability $$P(R_a=0,R_b=1|Y).$$ My idea is to sample from sample B with estimated selection probabilities $1-\hat{P}(R_a=1|Y)$. If I do this a large number of times, does the resulting sample have 'on average' the required selection probability? **Edit** : This feels very Bayesian. The question could maybe be put differently as: is it possible to simulate the bivariate distribution of two discrete variables from their marginal distributions?


title: Sample Balancing

text: Trying to figure out a way to balance small samples. Comparing ratios received from a survey, but sometimes I will have 25 survey responses and sometimes I will have 1. The ratios I look at might be the same but I think I need to take into consideration that one ratio is based on 25 surveys and the other is based on 1. Any way to do that with small samples? Thanks.


title: Is the difference between phantom and reference significant?

text: For 50 patients, 2 measurements were done: 1 on a phantom, 1 on a reference point. Ideally, the difference between the measurements would be zero. I want to test whether the difference is significant. Can anybody help me with this problem? Statistics isn't my strongest side. I hope I stated the situation well. Thanks.


title: How does the squared residual of a high leverage sample observation differ from the population's variance?

text: Is the high leverage sample observation considered a subset of the population, forcing its squared residual to be less than or equal to the population variance? Or am I misinterpreting the question? I'm trying to determine why the answer is false to the below question. ![Midterm Review Q](http://i.stack.imgur.com/XtsVk.png)


title: How different is training a factor graph with discriminative features?

text: Many of the people define graphical models with factors, each with 'conditional probability tables' (CPT) and perform inference on them. But more realistic case is when you can't define full probability tables, and instead you define functions of features with some probabilities (as in discriminative models): ![111](https://f.cloud.github.com/assets/2441454/2120733/c76de34a-91c1-11e3-969b-c846792a48fa.png) ![222](https://f.cloud.github.com/assets/2441454/2120738/115fd7e2-91c2-11e3-87a0-47a67fc9d834.png) Suppose we have an algorithm for inference on factor graphs defined with CPT. Implementation-wise, can I just use the same algorithm for inference on the above log-linear model?


title: How to compare the significance of two models from two different datasets?

text: I have two different regression models which I learned from two different data sets. Is there any statistical method which shows the significance of models based on the number of parameters and cross validation errors? The problem is that I used two different data sets and I worry that it might not make sense to compare two different models from two different data sets.


title: Constructing a pivot-based confidence interval

text: So I've been working on a problem in my probability class on which I have become stuck. It involves `X1 X2 ... Xn ~ Poisson(lambda)` 1 - We were instructed to show that `X_bar` was sufficient. Sol'n: I did this using factorization, and by finding `T(x) = Sum(x_i)/n`, which came from the exponent position above my lambda from the likelihood/pdf. It could've also been achieved via the exponential family method. 2 - We were then instructed to show that `X_bar^2` is approximately normal with mean `\lambda^2` and a variance depending on `\lambda` Sol'n: This was achieved using the central limit theorem and the delta method, which ended up giving me `N(0, 4*lambda^3)` **In step 3, however** , we are asked "to display an approximate pivot based on the result in part #2, and then produce a 95% confidence interval for `lambda^2`" We were given a hint, namely that               sqrt(n)(x_bar^2 - lambda^2) / sqrt (lambda * 4 lambda^2) ~ N(0,1)      I have, up until this point, no concept of what to do with a pivot or with this given info...I am looking for a number of degrees of freedom for a poisson table, which I think is wrong, but I've clearly got an issue with understanding what is being asked of me. If anyone could give me a nudge in the right direction, I'd be really grateful.


title: K-Means Clustering - Calculating Euclidean distances in a multiple variable dataset

text: I have just completed a simple exercise with 2 variables (X and Y) to understand how K-Means clustering works. The results look like this, ![Simple K-Means Clustering Result](http://i.stack.imgur.com/QAl9d.png) My questions is, if I have another column Z, how should the scatter plot be plotted to include the new variable Z? Would a 3 dimensional scatter plot chart be required? And how do I calculate the Euclidean distance between the coordinates for X,Y and Z? To clarify, I am not looking for a software like R or Weka to solve the problem but more on understanding the details and how the calculation works.


title: Good approximations for CDFs

text: As you probably know, the CDFs for many widely-used distributions are difficult, if not impossible, to express in closed form without the use of special functions   * The normal distribution uses $\mathrm{erf}$.   * Student's $t$-distribution requires the incomplete beta function or the Euler gamma function   * The $\chi^2$ distribution requires the Euler gamma function _and_ the lower incomplete gamma function. I could continue, but you get the idea. Unfortunately, all of these functions share the property of being impossible to compute. Does anybody know of a good (freely accessible) resource for approximating these various CDFs? Since I'm after p-values, I'm especially interested in approximations which are accurate for extreme values.


title: How to make an inference from one group to another

text: **Background:** I'm working on a statistics project that includes data for people who post a certain project online they need another to complete. Many people will apply to work on the project, but the poster will only choose one to complete it. There is a "back and forth" between the two parties until the project is completed. My job is to find all the projects that haven't been completed and to figure out whether they will ever be completed (ie the poster doesn't care about finishing the project, but it's still posted in the system - we want it gone). I've started by analyzing the data for "completed" projects to later compare to "not yet completed" projects. There are areas such as "how long to post a project after signing up", "how long to get a first response to a project", "how many applications until the poster makes a selection," etc. Here is a histogram of what each of the data columns looks like. Each seems to require a log transformation. (Ex: x is number of days to respond). ![Histogram](http://i.stack.imgur.com/7sRRq.png) **My main question is this:** Am I approaching this the correct way? Should I find the mean for and standard deviation for a variable such as "days to get a response" (from the "completed" data) and then if there's a value outside two standard deviations from the "not yet completed" data assume that project will not be completed? Or should I group all the variables together and do a multivariate analysis and compare the Mahalanobis distance of the "not yet completed" group? My problem with the multivariate approach, however, is that it's hard to think about the "ellipse around the data" (because it seems to be bounded by a square). Below is a plot variables against each other with their log transforms. ![LogTransform](http://i.stack.imgur.com/OJfmX.png) Or should I look at correlations? Any help with how to approach would be much appreciated. Thanks.


title: Good introduction book on small area estimation?

text: What would you guys suggest as a good introduction book to small area estimation?


title: -1 df for a saturated model with ML estimator?

text: I estimate a saturated / just-identified model (specifically: an APIM model, actor-partner-interdepence model) with the R package `lavaan`. This model is saturated and should have 0 $df$. This is true for the "GLS" generalized least squares estimator. With "ML" maximum likelihood estimation, however, it has **-1** $df$ according to the output. **Can somebody explain why there's a df of -1 with maximum likelihood estimation?**


title: Inconsistency in calculating the Calinski-Harabasz index for a given clustering in R

text: I am interested in determining the optimal number of clusters calculated by the PAM clustering algorithm using the Calinski-Harabasz (CH) index. To that end, I found 2 different R functions calculating CH values for a given clustering, but which returned different results: ?cluster.stats (in the fpc package), and ?index.G1 (in the clusterSim package). First one is called via:               pam.res <- pam(dist.matrix, 2, diss=TRUE)     ch1     <- cluster.stats(dist.matrix, pam.res$clustering, silhouette=TRUE)$ch      Second one is called via:               ch2 <- index.G1(t(dataframe), pam.res$clustering, d=dist.matrix)      Data may be found here: dataframe.RData, or here: dist.matrix.RData [dead links].   * **Can anybody explain the difference between these two CH index calculations to me?** Using `cluster.stats()`, the highest CH index is obtained for 2 clusters ($\approx32$); while using `index.G1()`, the highest CH index is obtained for 3 clusters ($\approx60$, and the value for 2 clusters is totally different from the previous, $\approx54$).   * **Which function is normally used to calculate the CH index?**


title: To chi square or not to chi square?

text: I am running a frequency count on the instances of powerless of language in the language of prosecuting attorneys and the same frequency count of the defendants language. I am asking if the frequency of the powerless language by women makes a favorable verdict and less powerless language by men creates a favorable verdict. I have 400 cases/transcripts. Below is my hypothesis: Is the use of powerless language used by defense attorneys in arguments independent of the verdict?               a. Is the use of powerless language used by female defense          attorneys in closing    arguments independent of the verdict?      Ha) The use of powerless language by a female defense attorney is          not independent of the verdict.     Ho) The use of powerless language by a female defense attorney is                        independent of the verdict        b. Is the use of powerless language used by male defense attorneys          in closing  arguments independent of the verdict?      Ha) The use of powerless language by a male defense attorney is          not independent of the verdict.     Ho) The use of powerless language by a male defense attorney is                          independent of the verdict.       Would I run two separate chi squares: one on all of the prosecuting attorneys and one on all the defending attorneys with the verdicts?


title: Restricted Boltzmann machines vs multilayer neural networks

text: I've been wanting to experiment with a neural network for a classification problem that I'm facing. I ran into papers that talk of RBMs. But from what I can understand, they are no different from having a multilayer neural network. Is this accurate? Moreover I work with R and am not seeing any canned packages for RBMs. I did run into literature that talks about deep learning networks which are basically stacked RBMs but not sure if it is worth the effort to implement them in R. Would anybody have any pointers? Thanks


title: How to model binary dependent data with temporal autocorrelation?

text: I am trying to model annual tree nut production using climate predictors. The nut data (dependent) is a binary timeseries (0,1 - representing unsuccessful and successful nut production), with one observation per year, and with 90 years of data and two missing years (88 onservations). The independent variables are monthly climate variables, including months in previous years (for example, Temp.July.t, Temp.July.t-1) I'm using R, and have an basic-intermediate knowledge of statistics. My problem is that the dependent data has strong temporal autocorrelation (nut production cannot be successful two years running). I'm looking for a pointer towards a technique that will allow me to deal with the autocorrelation in the binary data and create a statistical model that allows me to investigate the relationship between nut production and climate. Thank you.


title: Is the Moving Average of ARMA the same of Moving Average of Stock Market?

text: I'm studying time series prediction and I have some questions. Is the Moving Averages movel studied the methods of the ARMA family has the same concept as the methods studied in Moving Averages technical analysis stock trading? Is there any place where I can find implementations in `.Net` or `Java` of time series forecasting with Auto-Regressive `AR(p)`, Moving Average `MA(q)`, `ARMA(p, q)` and `ARIMA (p, d, q)`? Thank you


title: multivariate sorting / ranking

text: I have data for 1000 students' performance over 10 different tests on a scale of 0 -100 (a 1000 rows X 10 col matrix). I calculated the mean score and the associated coefficient of variation for each student. Now I wish to sort / rank students who consistently perform better, that is, have a high mean score and low coeff. of variation. Can anyone suggest a ranking scheme based on two objective variables (mean score and coeff. of variation). Thanks.


title: How do you report Kruskal Wallis one way analysis of variance and post hoc results following APA 6th standards?

text: How do you report Kruskal Wallis One Way Anova and post hoc results according to APA (American Psychological Association) 6th ed. standards? I have tried to find a template for reporting the results for the Kruskal-Wallis one way analysis of variance test and post hoc tests. I have found a number of descriptions but the templates all contained different types of information. Some included the full sample size, some included the sample size for the two groups being compared, some included the k-df, and some included the confidence intervals in text but I could not find anything consistent. The APA 6th edition’s examples are not specific to the Kruskal-Wallis and I don’t want to make assumptions. Each of the Kruskal Wallis analysis were significant at .000. I have a total of 392 post hoc comparisons: > 8 DAT groups = 28 comparison pairs between these groups.   >  15 SDRT tests. The SDRT is a Reading comprehension test. The items in this > test were divided into 14 sub-tests (plus total score). The sub-tests are > not all independent of each other (One item can be present in two or more > sub-tests). These represent different text conditions. Thus 28 group comparisons were run for all 15 SDRTs (total score and 14 sub- tests) which resulted in 392 analyses. I first reported the trends for the SDRT total score. Then I went through each SDRT sub-test to see if these trends changed when certain SDRT items are included in the assessment but not others. Thus the comparison is between the 8 DAT groups as well as between the 15 SDRTs. Could you suggest the best way to report this in text or refer me to some examples?


title: R: looking for "time" clusters in a data set

text: I am new to R and seeking some advice. I have a set (~20M) of data describing on which step a process did fail or succeeded:               name;star_time;step1;step2;step3     A;2011-07-22 11:02:44;passed;passed;failed     B;2011-03-06 19:33:06;passed;failed;failed     ...      I want to search for clusters of processes which have mostly succeed on stepX (for a specific month, ...) and visualize it (scatter plot). Then if possible I would like to find some patterns in data -> as process B did in 80% of the runs successfully ended in stepY after previous run of process Q that same day. Any advice where to start where to begin? I don't have an overly fast machine, so any efficiency advices are welcome. Thanks.


title: multicollinearity question, X and X change score

text: Consider survey data from surgeries. $Y$ represents observed surgical quality and is measured post-surgery; $X$ represents perceived surgical difficulty level and is measured pre and post surgery. It is desired to assess the relationship between $Y$ and $X_{pre}$ and also $Y$ on $X_{\delta} = X_{post} - X_{pre}$. However, since $X_\delta$ is derived from $X_{pre}$, we know that $X_{pre}$ and $X_{\delta}$ will be highly correlated. One option is to attempt to reduce such multicollinearity via centering.   Any thoughts on alternative strategies and pros/cons? This scenario sounds similar to the commonly discussed scenario of how to handle change scores via ANCOVA (e.g., Senn 2009, Stat Med) when we have change scores with respect to $Y$, but it is different since here we have no baseline $Y$ and a change score for $X$.


Is there a statistical measure for how much a variable fluctuates over time? For example a noise signal fluctuates a lot. However, if you would sort all values of the signal in time, you would have a signal that does not fluctuate at all. My first idea was to use the variance of the derivative, but I’m not sure whether this is the best measure. Another idea might be to use information entropy, but I’m not sure how to approach it that way.


title: Dirichlet density with just one x value and one alpha parameter

text: Does it make any sense to apply the Dirichlet density function to only one x value and therefore one alpha parameter? (this would be the result of some bin merging) I'm asking because it appears that gtool's implementation always returns 0 in such case, whereas this implementation doesn't. And this is the only case where their output they differs.               x=[0.94736842105263153];      alphas = [0.032034681077891182]     ======> implementation 1: 1.05372888806 ;      ======> gtool implementation: 0.      I'm very tempted to declare this result as non valid and re-merge bins differently.


title: What’s the most appropriate statistical analysis to compare predictions to observations?

text: I have a small set (say, 200 data points) of estimates about the time to complete a varied group of tasks, along with the actual times used to complete those tasks. These estimates are divided in four groups: 0:30, 1:00, 2:00, and 4:00. These groups indicate that the task should be completed in 30 minutes, 1 hour, and so on. The hypotheses that I’d want to test are:   * What is the correlation, if any, between the estimate and the difference with the actual times? In other words, are the smaller groups more precise and/or more reliable?   * Are the times significantly under- or overestimated?  PS _I don’t have a solid knowledge of statistics beyond the really really basic stuff. Maybe my question is meaningless, I don’t know. If that’s the case, sorry for the noise._


title: How to perform proper data mining on time-series data?

text: I have some daily data from city A, B, C. Values from city A are highly correlated with values from other cities for lag -1,-2,-3 and -4. I want to use Random Forest, SVM and ANN to predict values for city A. My idea is:   1. Split data into training and testing set.   2. Use the formula: valueA ~ valueB-1 + valueB-2 + valueB-3 + valueB-4 + valueC-1 + valueC-2 + valueC-3 + valueC-4   3. Try different methods (Random Forest, SVM and ANN) on training set and use createTimeSlices to cross-validation for model training and parameter tuning, like in this example - http://stackoverflow.com/a/22338029/2602477   4. Evaluate obtained model on the testing set (using R2, RMSE, etc.) My questions are:   * How to properly split data into training and testing set? I'm not sure, but boostraping or k-fold cross-validation doesn't sound right.   * Are given formula is correct for this case?


title: Large number of correlations with MRI region of interest (ROI) variables: What adjustment should I perform?

text: In our group, we have measured grey matter volumes of about 20 regions of interest (ROI) and frequency of cannabis use. Now, we would like find out whether cannabis use is associated with grey matter volumes of any of the measured brain regions. Furthermore, we would like find out whether these associations differ in subgroups of our sample. For each of the 20 ROIs, we have performed an ANCOVA with the ROI as dependent variable and cannabis use, group, sex, age, and total grey matter volume as independent variables. We are mainly interested in the effects of cannabis. The other independent variables in the models serve as control variables. To investigate group specific associatons between cannabis use and ROIs, we have specified interaction terms between cannabis use and group. What kind of adjustment for multiple comparisons would you use in such an analysis? I have already adjusted my p-values with the Benjamini-Hochberg procedure, but I'm concerned that it is too conservative because the ROI variables and consequently the associations between ROIs and cannabis use are not independent. I guess it would be more appropriate to analyze these associations in a single model such that the non-independency of the associations could be taken better into account. I was thinking about a mixed effects models in which the regression coefficient for cannabis use is allowed to randomly vary between ROIs. A big advantage of such an approach would be that one could adjust not only the standard errors of the regressions coefficients but also the parameter estimates themselves (due to shrinkage). Andrew Gelman has proposed such an approach in the following papers: http://pps.sagepub.com/content/4/3/310.short http://www.stat.columbia.edu/~gelman/research/unpublished/multiple2.pdf Do you think this would be feasible in my situation? How would such a model be specified in R?


title: Can odds ratios increase if you include more variables?

text: Suppose one wants to compute the odds ratio of a disease for a person in group $A$ versus that of a person in group $B$. Suppose we consider the following: (i) include age in the odds ratio calculation, (ii) include age and height in the odds ratio calculation. Is it possible for the odds ratio for (ii) to be much greater than the odds ratio of (i)?


title: Interpretation for simple slope analysis for curvilinear regression with interaction effects

text: When regressing income ($Y$) on age ($X$) moderated by gender ($Z$), I not only find significant effects for age ($X$), age squared ($X^2$), gender ($Z$), the interaction of age and gender ($XZ$), and the interaction of squared age and gender ($X^2Z$). Could anyone help me with how to interpret these results? Specifically, how can I calculate simple slopes for the interaction effects of age squared and gender in SPSS?


title: Using $\chi^2$ to compare two Markov transition matrices

text: I have a set of observed Markov sequences for which I have calculated first and second order transition matrices: $$M={P}(X_i=x_i|X_{i-1}=x_j)$$ and $$M^2=P(X_i=x_i|X_{i-1}=x_j,\,X_{i-2}=x_i )$$. I then produced a series of Markov chains based on the transition matrix $M$ and hence recalculated the empirical transition matrix $M^2$ based on this new data calling it $\tilde M^2$ and would like to compare it against the expected $M^2$. Is the $\chi^2$ test appropriate here? ie. $\displaystyle \sum \dfrac{M^2_{ij}-\tilde M^2_{ij}}{\tilde M^2_{ij}}$?


title: Approximating integrals using Monte Carlo simulation in R

text: How do I approximate the following integral using MC simulation? $$ \int_{-1}^{1} \int_{-1}^{1} |x-y| \,\mathrm{d}x \,\mathrm{d}y $$ Thanks! **Edit** (Some context): I am trying to learn how to use simulation to approximate integrals, and am getting some practice done when I ran into some difficulties. **Edit 2+3** : Somehow I got confused and thought I needed to split the integral into separate parts. So, I actually figured it out:               n <- 15000     x <- runif(n, min=-1, max=1)     y <- runif(n, min=-1, max=1)     mean(4*abs(x-y))


title: Eighth order moment

text: I read Nonlinear Dimensionality Reduction by Lee and Verleysen [Google Books] and came across the following theorem (p. 8): > Let $\mathbf{y}$ be a $D$-dimensional vector $[y_1, \ldots, y_d, > \ldots,y_D]'$; all components $y_d$ of the vector are i.i.d. with a finite > eighth order moment. Can anyone please explain what the authors mean by "finite eighth order moment"?


title: Interpretation of spatial Gcross plots

text: On page 192 of Analysing spatial point patterns in 'R' (Baddeley 2011), there are plots of the Gcross function for the amacrine dataset. I am looking for an interpretation of the plot. off/off and on/on are far below the CSR line, but on/off and off/on cross the CSR line - what does this mean please?


title: Comparing performance after CV over a single dataset : does it violate the independence assumption?

text: Cross validation lets us get multiple estimates of the generalization performance of an algorithm when the performance cannot be estimated for each example independently. (As with an F1-score, we need to estimate it over all of the data.) But when we use a t-test to determine whether or not there is a difference in performance between the two algorithms, we need to make an assumption that the estimates are independent. Do estimates derived from cross-validation violate this independence assumption?


title: co-twin study regression modelling

text: I want to conduct a regression analysis on a dataset of identical twins. There is a paper by Carlin et. al. 2005 which reviews the models that are designed for twin studies and includes the Stata codes in the appendix. The issue that I am facing is that I am not sure if i'm implementing them correctly. The variables that i'm studying are   * _jobClass_ as my dependent dichotomous variable, where \begin{equation} jobClass = \left\{ \begin{array}{l l} = 1 \ \text{ if whiteCollar} \\ = 0 \ \text{if blueCollar} \end{array} \right. \end{equation}   * _education_ as my categorical and ordered, independent variable (i.e. where 1 is the lowest education and 5 is the highest) One thing I did was that i also created a dichotomous independent variable which I called _uniEducation_. This variable is derived form the _education_ variable and is simply \begin{equation} uniEducation = \left\{ \begin{array}{l l} = 1 \ \text{ if person has university education} \\ = 0 \ \text{if otherwise} \end{array} \right. \end{equation} I want to conduct the study on each of the independent variables separately (first on _education_ , and then on uniEducation) I am interested in running first a regression without taking into account twin effects, (i.e. as in model 1D in the paper stated above). I then want to run a regression with the twin effects (i.e. as in model 2D in the paper) I am not really familiar with these types of regressions so I don't know exactly how to write the code in Stata. I know the codes are already provided in the appendix of the paper but i was told that I should run the regressions and include the fixed effects ",fe" at the end of the regression because twins have fixed effects. However I am not sure if the codes, as presented in the paper, already take that into account. I would greatly appreciate as detailed an explanation as possible, and how to exactly write the code in Stata. Thank you.


title: Moving from support vector machine to neural network (Back propagation)

text: I'm working with text recognition and currently I'm using support vector machine method. I would like to try with neural network also. I read a few documents about how neural network works, but the theory is quite heavy and I don't know exactly how will it apply to my case. So it would be good if someone can help me to make it clear, especially with the neural network's architecture.   * Currently, in SVM, I have 200 features (divided into 4 main categories), which is used to recognize the text. If I move to neural network, With 200 features, does it mean that I will have 200 neutrons in the input layer?   * With 200 features, how will that result in the architecture of the neural network (in term of numbers layers (hidden layer) and neutrons)?   * In SVM, I have one class classification (basically, true and false) and multi-class classification (labels), how this difference will apply to the output layer of the neural networks? And I also have a few general questions :   * What will help to decide the number of hidden layers and the number of neutrons inside each hidden layer?   * Does the number of hidden layer relate to the accuracy ? I'm new to neural network so it would be great if you can explain to me in an understable way. :) Thank you very much.


title: Repeated measures within and between conditions

text: I tested reactions to pictures under five conditions (neutral, high status, low status, high adhesion, low adhesion), using an EEG of 16 channels. The pictures were presented in random order for 15 times under each condition. So at best I have data from 75 trials for each subject in each channel. But as some electrodes sometimes made me some problems, I can't use the data from all of the EEG channels for all of the subjects. And as an EEG is sensitive to movements of the head, especially to eye movements, I furthermore had to reject some of the trials in all of the EEG channels for some of the subjects. So, on the whole, I neither have the same number of channels for each subject nor the same number of trials in the lasting channels for all of the subjects and therefore an unbalanced design. I split-up the data in 12 time segments per EEG channel and condition and I'm interested in the divergent effects of the conditions on the EEG waves in each of these time segments. So, I want to compare the data (trials) from all the subjects for condition "status high" in channel 1 for time segment 1 with the data from all these subjects for condition "status low" in channel 1 for time segment 1. And the data from all the subjects for conditions "status high" in channel 1 for time segment 2 with the data from all these subjects for condition "status low" in channel 1 for time segment 2. And so on, until the pairwise comparison for time segment 12 in that first channel. And then, I want to do the same with channel 2 and channel 3 and so on, until channel 16. And then, I do the same with "adhesion high" vs. "adhesion low" for all of the channels. Could I run these tests with the help of a Linear Mixed Model? What would be the fixed and random effects? Or would a multivariate approach be the better choice? I would be so grateful, if you could sacrifice a little bit of your time for my problem. Please excuse my maladroit wording.


title: What are RMSE SD and Rsquared SD metrics in resampling results using R package:caret?

text: I've been doing predictive modelling with R package `caret`. When resampling regression models, I get the traditional RMSE and Rsquared metrics, but also RMSE SD and Rsquared SD, for which I haven't found explanation in the manuals or documentations. Please, could anyone enlighten me? Reproducible code:               library(mlbench)     data(BostonHousing)     library(caret)          ctrl <- trainControl(method = "cv", number = 2)     lmFit <- train(medv ~ ., data = BostonHousing, method = "lm", trControl = ctrl)     print(lmFit)      outputs:               Linear Regression           506 samples     13 predictors          No pre-processing     Resampling: Cross-Validated (2 fold)           Summary of sample sizes: 254, 252           Resampling results          RMSE  Rsquared  RMSE SD  Rsquared SD     4.91  0.721     0.202    0.00304               Cheers Maria


title: How to get only desirable comparisons from post-hoc

text: I have some experimental data that I'm trying to analyze. I have 1 response variable and 3 explanatory variables (these are factor variables). The explanatory variables are the presence of a disease(positive and negative), a genetic profile (X and Y), and whether or not an MRI contrast agent was given (YES and NO). Structure of data looks like this:                    measurement   profile  disease contrast     1    -1.76269      X        NEG       YES     2    -0.34492      X        NEG       NO     3     0.57455      X        POS       YES     4     2.16539      X        POS       NO                 .      .          .        .                 .      .          .        .                 .      .          .        .     77   -1.76269      Y        NEG       YES     78   -0.34492      Y        NEG       NO     79    0.57455      Y        POS       YES     80    2.16539      Y        POS       NO      I looked into using ANOVA for this analysis but the post hoc **Tukey HSD looks at all possible combinations** of the explanatory variables so **it makes far more comparisons than I actually care about.** We have some specific hypotheses that, e.g.: X.NEG.NO will differ from Y.NEG.NO, X.NEG.NO will differ form X.NEG.YES, X.NEG.NO will differ from X.POS.NO, etc. _**(notice that each compared group "consist" from interaction of all three variables)_** How to get only some specific comparisons from TukeyHSD? Is appropriate to make Is there a better approach?               Reproducible example:          my.data <- data.frame(measurement = rnorm(80),                           my.profile = rep(c("X","Y"), each = 40),                           my.disease = rep(c("NEG","NEG","POS","POS"), times=20),                           my.contrast = rep(c("NO","YES"), times = 40))


title: Probability of getting an equal number of ones and twos from 10 dice rolls

text: If a die is rolled 10 times, how do you calculate the probability that an equal number of ones and twos occur? Let A be the number of ones, B the number of twos I first thought it might be equal to $\ (P(A=1) \times P(B=1)) + (P(A=2) \times P(B=2)) + ... + (P(A=10) \times P(B=10))$ but this can't work because $\ P(A=6) \times P(B=6) $ onwards is impossible. I know I could just stop after 5 each but I don't think that's the way to do it.


title: What is a full conditional probability?

text: I know what is a conditional probability distribution. But what is exactly _full conditional probability_?


title: How to get prediction for a specific variable in WinBUGS?

text: I am a new user of WinBUGS and have one question for your help. After running the following code, I got parameters of `beta0` through `beta4` (stats, density), but I don't know how to get the prediction of the last value of `h`, which I set to `NA` to model in the code. Does anyone can given me a hint? Any advice would be greatly appreciated. * * *               model {     for(i in 1: N) {     CF01[i] ~ dnorm(0, 20)     CF02[i]  ~ dnorm(0, 1)     h[i] ~ dpois (lambda [i])     log(lambda [i]) <- beta0 + beta1*CF03[i] + beta2*CF02[i] + beta3*CF01[i] + beta4*IND[i]     }     beta0 ~ dnorm(0.0, 1.0E-6)     beta1 ~ dnorm(0.0, 1.0E-6)     beta2 ~ dnorm(0.0, 1.0E-6)     beta3 ~ dnorm(0.0, 1.0E-6)     beta4  <- log(p)     p ~ dunif(lower, upper)     }          INITS     list(beta0 = 0, beta1 = 0, beta2 = 0, beta3 = 0, p = 0.9)          DATA(LIST)     list(N = 154, lower = 0.80, upper = 0.95,          h = c(1, 4, 1, 2, 1, 2, 1, 1, 1, 3, 3, 0, 0, 0, 2, 0, 1, 0, 4, 2,     3, 0, 2, 1, 1, 2, 2, 2, 3, 4, 2, 3, 1, 0, 1, 3, 3, 3, 1, 0, 1,     0, 5, 2, 1, 2, 1, 3, 3, 1, 1, 0, 2, 2, 0, 3, 0, 0, 3, 2, 2, 2,     1, 0, 3, 3, 1, 1, 1, 2, 1, 0, 1, 2, 1, 2, 0, 2, 1, 0, 0, 2, 5,     0, 2, 1, 0, 2, 1, 2, 2, 2, 0, 3, 2, 1, 3, 3, 3, 3, 0, 1, 3, 3,     3, 1, 0, 0, 1, 2, 1, 0, 1, 4, 1, 1, 1, 1, 2, 1, 3, 0, 0, 1, 1,     1, 1, 0, 2, 1, 0, 0, 1, 1, 5, 1, 1, 1, 3, 0, 1, 1, 1, 0, 2, 1,     0, 3, 3, 0, 0, 1, 2, 6, NA),          CF03 = c(-1.575, 0.170, -1.040, -0.010, -0.750,     0.665, -0.250, 0.145, -0.345, -1.915, -1.515,     0.215, -1.040, -0.035, 0.805, -0.860, -1.775,     1.725, -1.345, 1.055, -1.935, -0.160, -0.075,     -1.305, 1.175, 0.130, -1.025, -0.630, 0.065,     -0.665, 0.415, -0.660, -1.145, 0.165, 0.955,     -0.920, 0.250, -0.365, 0.750, 0.045, -2.760,     -0.520, -0.095, 0.700, 0.155, -0.580, -0.970,     -0.685, -0.640, -0.900, -0.250, -1.355, -1.330,     0.440, -1.505, -1.715, -0.330, 1.375, -1.135,     -1.285, 0.605, 0.360, 0.705, 1.380, -2.385, -1.875,     -0.390, 0.770, 1.605, -0.430, -1.120, 1.575, 0.440,     -1.320, -0.540, -1.490, -1.815, -2.395, 0.305,     0.735, -0.790, -1.070, -1.085, -0.540, -0.935,     -0.790, 1.400, 0.310, -1.150, -0.725, -0.150,     -0.640, 2.040, -1.180, -0.235, -0.070, -0.500,     -0.750, -1.450, -0.235, -1.635, -0.460, -1.855,     -0.925, 0.075, 2.900, -0.820, -0.170, -0.355,     -0.170, 0.595, 0.655, 0.070, 0.330, 0.395, 1.165,     0.750, -0.275, -0.700, 0.880, -0.970, 1.155, 0.600,     -0.075, -1.120, 1.480, -1.255, 0.255, 0.725,     -1.230, -0.760, -0.380, -0.015, -1.005, -1.605,     0.435, -0.695, -1.995, 0.315, -0.385, -0.175,     -0.470, -1.215, 0.780, -1.860, -0.035, -2.700,     -1.055, 1.210, 0.600, -0.710, 0.425, 0.155, -0.525,     -0.565),          CF02 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA,     NA, NA, NA, NA, NA, NA, 0.38, 0.06, -0.94,     -0.02, -0.28, -0.78, -0.95, 2.33, 1.43, 1.24, 1.26,     -0.75, -1.5, -2.09, 1.01, -0.05, 2.48, 2.48, 0.46,     0.46, -0.2, -1.11, 0.52, -0.37, 0.58, 0.86, 0.59,     -0.12, -1.33, 1.4, -1.84, -1.4, -0.76, -0.23,     -1.78, -1.43, 1.2, 0.32, 1.87, 0.43, -1.71, -0.54,     -1.25, -1.01, -1.98, 0.52, -1.07, -0.44, -0.24,     -1.31, -2.14, -0.43, 2.47, -0.09, -1.32, -0.3,     -0.99, 1.1, 0.41, 1.01, -0.19, 0.45, -0.07, -1.41,     0.87, 0.68, 1.61, 0.36, -1.06, -0.44, -0.16, 0.72,     -0.69, -0.94, 0.11, 1.25, 0.33, -0.05, 0.87, -0.37,     -0.2, -2.22, 0.26, -0.53, -1.59, 0.04, 0.16, -2.66,     -0.21, -0.92, 0.25, -1.36, -1.62, 0.61, -0.2, 0,     1.14, 0.27, -0.64, 2.29, -0.56, -0.59, 0.44, -0.05,     0.56, 0.71, 0.32, -0.38, 0.01, -1.62, 1.74, 0.27, 0.97,     1.22, -0.21, -0.05, 1.15, 1.49, -0.15, 0.05, -0.87,     -0.3, -0.08, 0.5, 0.84, -1.67, 0.69, 0.47, 0.44,     -1.35, -0.24, -1.5, -1.32, -0.08, 0.76, -0.57,     -0.84, -1.11, 1.94, -0.68),          CF01 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA,     NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,     NA, -0.117, -0.211, -0.333, -0.229, -0.272,     -0.243, -0.148, 0.191, -0.263, -0.239, -0.168,     -0.381, -0.512, -0.338, -0.296, 0.067, 0.104,     -0.254, -0.167, -0.526, -0.096, -0.43, 0.013,     -0.438, -0.297, -0.131, -0.098, -0.046, -0.063,     -0.194, -0.155, -0.645, -0.603, -0.374, -0.214,     -0.165, -0.509, -0.171, -0.442, -0.468, -0.289,     -0.427, -0.519, -0.454, 0.046, -0.275, -0.401,     -0.542, -0.488, -0.52, -0.018, -0.551, -0.444,     -0.254, -0.286, 0.048, -0.03, -0.015, -0.219,     -0.029, 0.059, 0.007, 0.157, 0.141, -0.035, 0.136,     0.526, 0.113, 0.22, -0.022, -0.173, 0.021, -0.027,     0.261, 0.082, -0.266, -0.284, -0.097, 0.097, -0.06,     0.397, 0.315, 0.302, -0.026, 0.268, -0.111, 0.084,     0.14, -0.073, 0.287, 0.061, 0.035, -0.022, -0.091,     -0.22, -0.021, -0.17, -0.184, 0.121, -0.192,     -0.24, -0.283, -0.003, -0.45, -0.138, -0.143,     0.017, -0.245, 0.003, 0.108, 0.015, -0.219, 0.09,     -0.22, -0.004, -0.178, 0.396, 0.204, 0.342, 0.079,     -0.034, -0.122, -0.24, -0.125, 0.382, 0.072, 0.294,     0.577, 0.4, 0.213, 0.359, 0.074, 0.388, 0.253, 0.167),          IND = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,     1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,     1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,     0, 0, 0, 0, 0, 0, 0, 0, 0, 0))


title: I have my data set, now what?

text: I have a basic understanding of basic statistics, but I believe I've gotten myself out of my depth. I have a data set with a dependent variable (time span) and three quantitative independent variables. There is also a qualitative independent variable (a type flag) but I think I can quantify it if I must. The quantitative independent variables all seem to be inversely proportional to the dependent variable. I want to see if I can determine which of these independent variables has the most influence over the dependent. I ran all data through regression tests using `statsmodels` for Python (ordinary least squares) but I get warnings indicating multicollinearity and $R^2 = 0.002$. I've also tried univariate linear regression of each independent but am not getting something that looks usable. Finally, I also have this problem of the qualitative independent. Again, I'd really like to know how these independents work together to influence the dependent, and I'd also like to know the degrees to which they influence but I'm clearly lost as to the methodology. **EDIT** > Warnings: > > 1 The condition number is large, 5.65e+06. This might indicate that there > are strong multicollinearity or other numerical problems. Screenshot of full output from `.fit().summary()` ![summary](http://i.stack.imgur.com/iwZvn.png) Data Profile: ![enter image description here](http://i.stack.imgur.com/mgYdB.png) I should also note that two of these independent variables (`maxtransfersize` and `buffercount`) repeat for every combination of `stripes` and `diskconfig` (<- this is the qualitative variable, and I think I have a way of quantifying that now). `backup_time_ms` is the dependent variable.


title: Is there a software package to fit grouped (aggregated) logistic panel regression models?

text: I have a dataset that looks like the following               id | time | trial  | success | other covariates      1    0  10    4         1    1  23    10          1    2  100   29          2    0  5     1         etc this is basically a panel aggregated logistic regression. Can someone point me towards a package I can use for grouped logistic PANEL regression? I am aware of PGLM package in R I tried to use stata with xtlogit but it wouldn't take grouped logistic with fweight options. What I am looking for is something like blogit in stata but something that can take panel data. I'd presume the function would be called xtblogit but it doesn't exist. Thank you!


title: Getting past independence assumptions in modeling the sum of random variables (application in education)

text: I'm trying to model a student's semester GPA $G$ as a random variable. Semester GPA is a weighted (based on credit hours) sum of a student's grade points (e.g. 0=F,1=D,2=C,3=B,4=A). We can consider the students grade points in each course $G_1,...,G_n$ as random variables, each with a particular distribution as well. It could be discrete, as it actually is, or it could be approximated by a continuous distribution (possible the beta distribution, since this is limited to an interval). If we assume independence of $G_1,...G_n$ then getting the distribution of $G$ is not so hard. But this is not so realistic. A few problems with this: students taking a lot of hours are likely to do worse overall than they would if they were just taking 1 of the courses, just because they have a limited amount of time. Additionally, students who do well in a particular course are likely to do well in a similar course. How would you approach incorporating this sort of dependence?


title: R package for combining factor levels for datamining?

text: All, Wondering if anyone has run across a package/function in R that will combine levels of a factor whose proportion of all the levels in a factor is less than some threshold? Specifically, one of the first steps in data preparation I conduct is to collapse sparse levels of factors together (say into a level called 'Other') that do not constitute at least, say, 2% of the total. This is done unsupervised and is done when the objective is to model some activity in marketing (not fraud detection, where those very small occurrences could be extremely important). I am looking for a function that will collapse levels until some threshold proportion is met. Does anyone know of anything? Thanks! Brian **UPDATE:** Thanks to these great suggestions I wrote a function pretty easily. I did realize though that it was possible to collapse levels with proportion < the minimum and still have that recoded level be < the minimum, requiring the addition of the lowest level with proportion > the minimum. Likely can be more efficient but it appears to work. The next enhancement would be to figure out how to capture the "rules" for applying the collapse logic to new data (a validation set or future data).               collapseFactors<- function(tableName,minPercent=5,fillIn ="RECODED" )     {         for (i in 1:ncol(tableName))             {                         if(is.factor(tableName[,i]) == TRUE) #process just factors                 {                               sortedTable<-sort(prop.table(table(tableName[,i])))                     numberToCollapse<-length(sortedTable[sortedTable<(minPercent/100)])                          if (sum(sortedTable[1:numberToCollapse])<(minPercent/100))                         {                             numberToCollapse=numberToCollapse+1 #add next level if < minPercent                         }                          if(numberToCollapse>1) #if not >1 then nothing to collapse                     {                         lf <- names(sortedTable[1:numberToCollapse])                         levels(tableName[,i])[levels(tableName[,i]) %in% lf] <- fillIn                     }                 }#end if a factor                       }#end for loop              return(tableName)          }#end function


title: Sample size for proportions in repeated measures

text: I'm trying to help a scientist design a study for the occurrence of salmonella microbes. He would like to compare an experimental antimicrobial formulation against a chlorine (bleach) at poultry farms. Because background rates of salmonella differ over time, he plans to measure % poultry w/salmonella before treatment, and after treatment. So the measurement will be the difference of before/after % salmonella for the experimental vs. chlorine formulas. Can anyone advise on how to estimate the sample sizes necessary? Let's say the background rate is 50%; after bleach it's 20%; and we want to detect whether the experimental formulation changes the rate by +/- 10%. thank you EDIT: What I'm struggling with is how to incorporate the background rates. Let's call them p3 and p4, the "before" salmonella rates for bleach and experimental samples, respectively. So the statistic to be estimated is the difference of differences: Experimental(After-Before) - Bleach(After-Before) = (p0-p2) - (p3-p1). To fully account for the sampling variation of "before" rates p2 and p3 in the sample-size calculation --- is it as simple as using p0(1-p0)+p1(1-p1)+p2(1-p2)+p3(1-p3) wherever there's a variation term in the sample-size equation? Let all samples sizes be equal, n1 = n2 = n.


title: Incremental learning for LOESS time series model

text: I am currently working on some time series data, I know I can use LOESS/ARIMA model. The data is written to a vector whose length is 1000, which is a queue, updating every 15 minutes, Thus the old data will pop out while the new data push in the vector. I can rerun the whole model on a scheduler, e.g. retrain the model every 15 minutes, that is, Use the whole 1000 value to train the LOESS model, However it is very inefficient, as every time only one value is insert while another 999 vlaues still same as last time. So how can I achieve better performance? Many thanks


title: Goldfarb Idnani quadratic solver

text: I am implementing the Support Vector Regression (SVR) algorithm by means of quadratic programming. In order to do that, I am using an optimization library that contains a quadratic solver based on the Goldfarb Idnani active set quadratic optimization method. However, the library seems not to be working the way I was expecting. Once the optimization process is finished, not all the constraints are active, so there are some variables for which the value is out of the valid range of values in an SVR problem. I implemented SVR in the past in Matlab by means of the quadprog method, and I did not have such issues. Of course, I reviewed and debugged my own code. My question is: is it normal that, after finishing the Goldfarb Idnani algorithm, some of the constraints are not active? Or should I assume that the library is not working properly? PS: I have read the original paper[1], and from it I guess that all constraints should be active at the end of the process. [1] "A numerically stable dual method for solving strictly convex quadratic programs", D. Goldfarb, A. Idnani, Mathematical Programming 27, pp. 1–33, 1983. Available at http://www.javaquant.net/papers/goldfarbidnani.pdf.


title: Missing at Random Data in GEE

text: For a continuous outcome being analyzed using GEE with a linear link, you have assurance that standard errors and point estimates are consistent with a first order trend regardless of distribution of outcome, heteroscedasticity, and mild non-linearity problems. Point estimates from the GEE are the same as those obtained from maximum likelihood (OLS), but the standard error estimates are the HC sandwich based errors and thus swamp up mild bits of classical model assumption violations. In longitudinal analyses where attrition depends upon measured variables (e.g. age), you know that the so-called "missing data mechanism" is missing at random (not missing COMPLETELY at random, per Little, Rubin 2002) and, further, that maximum likelihood estimates "are not biased" due to the factorization of the likelihood including the missing data indicator and unobserved likelihood contribution due to measured rows. My questions are:   1. For ML estimates, are complete case analyses considered _efficient_?   2. For GEE with linear link, are estimates somehow _biased_ even though they're the same as those obtained from ML?    3. Is the real problem that SEs from GEE with linear link are not guaranteed to be consistent? More so than is attributable to effective sample size loss due to complete case analysis?   4. Does weighting promise to help remedy the the SEs above and beyond effective sample size loss due to complete case analysis if there are other reasons why the GEE would be "wrong" in this case?


title: is it possible to check Czuprow's coefficient of the convergence in R package?

text: just wondering - is it possible to calculate Czuprow's coefficient of the convergence in R package? if so, is there some dedicated package for that purpose? thanks in advance


title: Does "correlation" also mean the slope in regression analysis?

text: I'm reading a paper and the author wrote: > **The effect of A,B, C on Y was studied through the use of multiple > regression analysis. A,B,C were entered into the regression equation with Y > as the dependent variable. The analysis of variance is presented in Table 3. > The effect of B on Y was significant, with B correlating .27 with Y.** English is not my mother tongue and I got really confused here. First, he said he would run a regression analysis, then he showed us the analysis of variance. Why? And then he wrote about the correlation coefficient, is that not from correlation analysis? Or this word could also be used to describe regression slope?


title: construct/load dataset that performs better with diffusion kernel than other kernel

text: I'm looking for a dataset on which a diffusion kernel (also called heat kernel), used via SVM, would get better accuracy than other kernels for the classification task. I want to use such a dataset to make sure my implementation of the diffusion kernel is working, and to get a better understanding of the diffusion kernel (does it classify as "similar" graphs that have the same expression on the same zones (doesn't care about the level expression as long as it's approximatively the same value on some part of the graph) ? Or can it understand substructure of the graphs, and say that for instance if 2 out of 3 connected features is present, it belongs in the same class ?). Does anyone know of such a dataset, or know how to construct it ? Thanks! Jess PS : For now, I've tried creating my own dataset in R, but the linear and gaussian kernel get similar results for classification. Here was the idea behind it : I'm trying to predict whether a pokemon is dragon type, based on this underlying graph : ![feature graph](http://i.stack.imgur.com/zTP2H.jpg) Here is the dataset, in R :                   charizard = c(70, 100, 80, 100, 80, 10, 95)         dragonite = c(80, 80, 80, 60, 20, 20, 50)         tyranitar = c(90, 0, 100, 80, 40, 5, 70)         altaria = c(30, 95, 10, 80, 60, 70, 80)         aerodactyl = c(65, 100, 100, 50, 60, 10, 40)         salamence = c(100, 100, 100, 100, 90, 0, 70)         rayquaza = c(100, 0, 100, 100, 80, 10, 40)         lugia = c(100, 80, 90, 80, 70, 40, 30)         gyarados = c(80, 0, 90, 80, 70, 2, 75)         hydreigon = c(95, 90, 90, 90, 70, 15, 70)         dratini = c(20, 0, 90, 60, 20, 60, 70)         hooh = c(100, 90, 5, 90, 70, 20, 35)              vibrava = c(60, 90, 20, 10, 100, 20, 50)         entei = c(95, 0, 5, 100, 60, 10, 60)         scyther = c(40, 40, 60, 0, 100, 20, 75)         magikarp = c(0, 0, 40, 0, 20, 20, 10)         butterfree = c(30, 100, 10, 5, 40, 90, 90)         milotic = c(90, 0, 100, 20, 80, 80, 90)         ponyta = c(15, 0, 30, 70, 70, 80, 60)         blastoise = c(90, 0, 40, 10, 40, 50, 80)         alakazam = c(95, 0, 60, 40, 100, 20, 60)         sudowoodo = c(50, 0, 60, 5, 30, 40, 40)                   train = rbind(charizard, dragonite, tyranitar, altaria, aerodactyl, salamence, rayquaza,                       vibrava, entei, scyther, magikarp, butterfree, milotic, ponyta)              test = rbind(lugia, gyarados, hydreigon, dratini, hooh,                      blastoise, alakazam, sudowoodo)              cl_train = as.factor(c(rep('dragon', 7), rep('not_dragon', 7)))              cl_test = as.factor(c(rep('dragon', 5), rep('not_dragon', 3 )))      Please note that I'm fully aware that Ho-oh is not a dragon pokemon, and neither is tyranitar nor lugia. It's just an example (and they kind of look like a dragon...). As mentionned earlier, this small dataset gets perfect accuracy for linear kernel, for instance, and not for the diffusion kernel, so it's not what I'm looking for even though I tried to buid it for this purpose...


title: Can you calculate a AIC value using the non-linear maximization (nlm) minimum value in R?

text: So the formula for AIC is: AIC = 2k - 2ln(L) L is the maximized value of the likelihood function. I'm modeling oxygen data in R using Non-Linear Minimization (nlm) of a maximum likelihood estimation (mle). By using the minimum command it produces the value of the estimated minimum of the function. This seems to be a rough measure of how well my model fits the data. I'm hoping to use this value as L so I can calculate AIC. Is this even possible? I don't understand very well what the nlm function is actually doing so as to know how well it fits with the likelihood function. Is there a better way to get the maximized value of the likelihood function from nlm? thanks


title: Practical Implementation of Gibbs Sampling in Latent Diriclet Allocation

text: In the collapsed Gibbs sampling version of LDA, the posterior distribution of topic assignments for each word is sampled. From what I have read (e.g. http://people.cs.umass.edu/~wallach/courses/s11/cmpsci791ss/readings/griffiths02gibbs.pdf) it looks like people just implement the Gibbs sampling procedure by letting it run until we are sure that we are sampling from the posterior then take the last sample as a good topic assignment. But shouldn't we be a bit more rigorous than this? For example, what if the sample we take happens to be from a part of the posterior with low probability (this would be unlikely but possible)? Wouldn't it be better to do something like take a Bayes estimator (or average, or mode) of the samples to make sure that we are using a topic assignment that is actually likely?


title: Covariates in Path model

text: I am building a path model which links different psychological constructs. This is the replication of a model I have previously tested in another sample. The model fits the data well. When however I try to add an additional variable as covariate on the DV (i.e., a variable with an arrow pointed to the DV) the model does not fit the data anymore (all indices decrease, chi-squared becomes significant). What does that mean conceptually? Is the model wrong, because adding the covariates make it not fitting anymore? Or I can justify not using this covariate just because once added the model does not fit the data anymore? Using simple multiple regression, my focal IV remains significant also after adding the covariate (which also remains significant). I am aware this is a very easy question - but I am new to statistics and these approaches! thank you for your answers! Additional Information: ![This is the path diagram](http://i.stack.imgur.com/YYcWL.png) T2R and T2N predict T2m which in turn predict T2c and T2C (the mediators), which in turn predict T2P. T2S is the covariate. When the covariate is added to the model, the model does not fit the data anymore


title: Time series pattern identification - SVD/SSA?

text: I've looked over other posts regarding time series data, and am unsure if the mentioned methods would apply to what I'm trying to do, since I'm not familiar with pattern analysis methods: I have time series data from two conditions. I want to identify a pattern in each condition, and then write a piece of code (in R, Matlab, etc.) that will alert me when one of these patterns re-appear in a time series (in real-time). What avenue of time series analysis should I pursue?


title: Generate two variables with precise pre-specified correlation

text: I want to generate two variables with (pseudo-) random numbers with an exact pearson's r. How do I do that? Python and/or R solutions would be nice! I am able to generate random data that approximates a pre-specified r in python in the following way. I'm not searching for approximations but for data that with an exact pre-specified r, i.e. with r=0.83000 in the example below:               samples = 200     r = 0.83          # Generate pearson correlated data with approximately cor(X, Y) = r     import numpy as np     data = np.random.multivariate_normal([0, 0], [[1, r], [r, 1]], size=samples)     X, Y = data[:,0], data[:,1]          # That's it! Now let's take a look at the actual correlation:     import scipy.stats as stats     print 'r=', stats.pearsonr(X, Y)[0]      The motivation for knowing r is that I'm testing out (bayesian) statistical models that can infer r from data and they are a lot easier to evaluate when r is well specified. **SOLUTION:** thanks to Greg Snow for pointing out the empirical=TRUE command in mvrnorm (multivariate random normal stuff)! Here's the explicit code:               samples = 200     r = 0.83          library('MASS')     data = mvrnorm(n=samples, mu=c(0, 0), Sigma=matrix(c(1, r, r, 1), nrow=2), empirical=TRUE)     X = data[, 1]  # standard normal (mu=0, sd=1)     Y = data[, 2]  # standard normal (mu=0, sd=1)          cor(X, Y)  # yay!     cor(X*0.01 + 42, Y*3 - 1)  # Linear transformations of X and Y won't change r.          Y = Y*3 + 1


title: use of t-test to compare performance of algorithms

text: I need a little bit guidance. I have to compare the classification performance of multiple algorithms using simple or paired t-test. Let's say I have four datasets (A,B,C) with training and test samples. I am running 3 algorithms (SIFT, SURF, ORB) and compute the classification accuracy such 0.9 means 90% of images correctly match from the test dataset. Let' say I get the following table: Dataset (A,B,C,D)   * SIFT (0.90, 0.84, 0.90,0.45)   * SURF (0.84, 0.67, 0.45,0.34)   * ORB (0.34,0.45,0.45,0.23) Can you please guide me how I can compare the performance of these algorithms using some statistical analysis such as simple t-test? Any guidance will be really appreciated. Thanks.


title: 2D multivariate normal coverage probability

text: I have a moving platform which moves in orthogonal axis to move a light detector to take account of positional uncertainty. I have a systematic error in space that is well described by a circular normal distribution (actually its not quite circular but the covariance matrix is almost exacty diagonal). I can think of two ways I could determine the amount by which my platform needs to move. I can either look at the probability as a function of radius so that my platform needs to move by 2.15 times the standard deviation to achieve 90%, this was my first though. The motors on the platform however move in orthogonal axes so I also though that I could look at the marginal distributions in which case for each axis the platform needs to be able to move by 1.64 times the standard deviation to achieve 90% coverage. As it happens this isn't really a concern as the range is much greater than required but it got me thinking. If I had a positional uncertainty distribution that was a circular normal distribution which method should I use to detemine the amount of movement I need to cover 90% of events. If I think just about the motion in a single axis then, thinking radially, I need to move by 2.15 times the standard deviation however thinking 'marginally' then I only need 1.64 times the standard deviation. Which method is correct?


title: Error computing Robust Standard errors in Panel regression model (plm,R)

text: I am using the plm library to run fixed effect regressions and the sandwich,lmtest libraries to compute robust standard errors. I have no problem running the regressions, but in some instances when I go to compute the standard errors I get the following error:               library(plm)     library(sandwich)     library(lmtest)          fe_reg <- plm(y ~ x + I(x^2)+factor(date), data=reg_data, index=c("id","date"), model="within")     coeftest(fe_reg, vcov.=vcovHC(fe_reg, type="HC1"))          RRuntimeError: Error in solve.default(crossprod(demX))       system is computationally singular: reciprocal condition number = 1.84726e-25      I do not have any problems computing the coefficients or the "normal" standard errors (ie homoscedastic). Also, I have no problem computing the robust standard errors when I omit the quadratic term:               fe_reg <- plm(y ~ x +factor(date), data=reg_data, index=c("id","date"), model="within")        Anyone have any idea what is going on? If the design matrix were singular then the coefficients should not have been computed, so I don't understand where the issue is coming from when computing the standard errors. Thanks!


title: Best way to compare distributions with spikes

text: Following this question, I am interested on the best way to display four distributions with spikes. Indeed, I have four distributions with circa 75% of elements in it at the maximum value (1). I am searching for a method that is able to display graphically which distribution is more clustered near the maximum value. I thought about boxplots, but, as you can see in the plot below, it is not easy to decide which one of the second and the fourth boxes has the most quantity of data clustered near the maximum value (to the right). ![boxplot](http://i.stack.imgur.com/cgkz8.png) I thought about violin plots, but, with data with such spikes I cannot understand how the density of the violins are computed, because at least 75% of data has the maximum value. ![violin](http://i.stack.imgur.com/AaZaE.png) I tried to use a barplot. But even if you use fine binning near the maximum value, you have always the 75% at top and it is difficult to me to compare these four distributions. ![barplot](http://i.stack.imgur.com/legoI.png) I thought about qqplot between each couple of distributions, but I am unsure if this method can be used to the scope of defining which distribution is better clustered to the maximum. ![qqplot](http://i.stack.imgur.com/kNQWa.png) Is there any other data visualization method to represent this data in such a way that I can decide which of the four distribution is more clustered near the maximum value? I don't want to use the mean because can be affected by a far outlier. Neither the median because here the median is always the maximum value. I was thinking about counting how many points are above a certain threshold, but, how to define a good threshold? Thanks


title: Clustering using longitude, latitude, and some other variable

text: I am hoping someone can point me in the right direction with this problem I am having. I am trying to cluster geographical areas (basically using latitude and longitude as the zip code centroid) which have the highest number of deaths in a given state. But, the added constraint is that those clustered areas have to be touching. For example, if California has 1000 deaths and Nevada has 1000 deaths, and all other states have 0, then California and Nevada should be cluster 1 while all other states should be cluster 2. However, if California has 1000 deaths and Texas has 1000 deaths, and all other states have 0, then California should be cluster 1, Texas should be cluster 2, while all other states should be cluster 3. This is a bit of a simplification because I'm trying to do this exercise at the zip code level. However, the point is that I want proximity to matter. Does anyone know of a method that can achieve this result?


title: Paired difference test

text: I am predicting missing values by using two different methods (e.g., data imputation methods). Following I am interested to test if the predicted values, obtained by means of two different methods, are significantly different. The missing values are not normally distributed and lay in the range between 0 and 1. What method should I use? Thank you for your help.


title: Validating residual plot count data (different levels)

text: I am studying the distribution of a marine species using the number of sightings as a dependent variable. When I am trying to validate the plots of the best model I am getting a non-usual pattern, and then I decided to get the null model with different distributions, and I don’t know what is causing the different levels on the residuals vs fitted values. ![enter image description here](http://i.stack.imgur.com/Q3xCK.jpg)               null<-glm.nb(Sightings~offset(logEffortScale),data=l, link=log)      ![enter image description here](http://i.stack.imgur.com/eoVWG.jpg)               null<-glm.nb(Sightings~offset(logEffortScale),data=l, link=log)          plot(residuals(null)~predict(null))      ![enter image description here](http://i.stack.imgur.com/a4ws9.jpg) **Edit: correct residuals vs predicted plot**               plot(residuals(null)~predict(null,type="response"))      ![enter image description here](http://i.stack.imgur.com/Ttqde.jpg)               m1 <- zeroinfl(Sightings ~ offset(logEffortScale), data = l)      ![enter image description here](http://i.stack.imgur.com/6KsFo.jpg)               m1b <- zeroinfl(Sightings ~ offset(logEffortScale), data = l, dist = "negbin")      ![enter image description here](http://i.stack.imgur.com/UOaEl.jpg)               m2 <- zeroinfl(Sightings ~ offset(logEffortScale), data = l, dist='poisson', link='logit')      ![enter image description here](http://i.stack.imgur.com/jxWrB.jpg)


title: Methods for probability estimates

text: I am trying to learn about methods for learning probability estimates. I understand how to use logistic regression for that, but I wonder if there are more complex nonlinear methods of the vein of neural networks and random forests that can give me similar results? Any help would be much appreciated.


title: How to perform exponential regression with more than one predictors using R

text: Fitting exponential regression of form $y=ab^{x_1}c^{x_2}$ Here $x_1$ and $x_2$ are predictors and $y$ is dependent variable how to calculate $a,b,c$ using R tool


title: Looking to generate fun graphs/data for basketball matches

text: I'm looking to generate as much interesting data as I can from some 2on2 basketball games. The only data I have recorded is the players on each team, and the final score. What interesting things can I learn from this data? I'm interested in what I can learn about the individual players, as well as the teams and player combinations. Can anyone provide a good starting point.


title: Probability in logistic regression

text: I come across a practice where the logit (Intercept + beta* variables) is once again input as a independant variable in a single variable logistic regression against dependant variable to get intercept ($I$) and a beta ($B_0$). Finally, probability is calculated: $\frac{1}{1+e^{-(I+B_0\text{logit})}}$. Why is this done -- since probability can be directly calculated from the logit in the first step?


title: How to group-center / standardize variables in R?

text: Functions I'm familiar with include scale from base R, rescale from ARM. Perhaps the best way would be to use some variant of apply, specifying one or more variables to use as grouping variables.


title: Deep and sparse vs. Thin and dense data collection for curve fitting?

text: I'm designing a toy experiment in which I'd like to model the 'attractiveness' of a credit card as a function of the interest rate it comes with:   * Attractiveness will be reported on a 10 point scale   * I can vary interest rate between 0% and 25%.    * In total, I have budget that will allow me to collect 1,000 data points, each of which will be a single human subject's attractiveness rating of the credit card for a single interest rate. My question: Is it better to test more interest rate levels over this range, with fewer observations at each interest rate, or more interest rate levels over the range, and collect fewer observations at each level? I think if I want to fit a straight line (e.g. using linear regression), then as long as I collect the same number of observations at each interest level, and 2 or more levels, it doesn't matter. However, I'm not sure what family of curves I'd like to fit, or where I'd find math to support this. Thanks!


title: Does the order of explanatory variables in logistic regression change the results?

text: I am curious about the consequences of changing the order of the explanatory variables in a binary logistic regression. In a recent series of logistic regressions I ran in SPSS, I found that changing the order of the explanatory variables ($z, y, x$ instead of $x, y, z$) resulted in different coefficient values and significance levels. Investigating further, I got the same results in `R` using the same orders. Clearly, shifting the predictors around matters in terms of the results—my question is how? (And yes, I checked to make sure I wasn't entering the variables stepwise.)


title: Reliability of mean of standard deviations

text: I've a question which probably is going to show my ignorance about statistics :). I have a large set of machines that produce iron bars of certain lengths. For each machine, I have ran experiments and have a list of lengths. From those I can calculate a mean and sample standard deviation. I don't really care about their means and I am mainly focused on the variation. Therefore, I basically only record a sample standard deviation per machine. I think the results of each machine follow a normal distribution. So far so good :) I now want to combine these variations into a single number. Therefore, I calculate the quadratic average of each machine variation, let's call it X. In the next step, I also would like to give an estimate for the spread around X. What is this number called and what's the best way to calculate it? I'm not sure it's related to the confidence interval of a standard deviation and I don't know whether the measurements are independent (a design fault would show up in all, a construction maybe only in some). **Example.** I'll try to clarify with an example. Suppose I measure 3 machines and find that they produce lengths of   M1: 100 +/- 7   M2: 120 +/- 8   M3: 130 +/- 9   where the numbers behind the +/-'s are the sample standard deviations of observed values on that single machine. As said before, I don't care about the means but only in the spread, so I define {X_1, X_2, X_3} = {7,8,9}. Their quadratic average is X = RMS(X_i) = $\sqrt{194}$ and I think of X as an indication of the average spread of a machine in my park. Suppose that I would have found {X_1, X_2, X_3} = {3,8,11}. Their quadratic average is the same $\sqrt{194}$, but the spread around it is obviously bigger. My confidence in the correctness of $\sqrt{194}$ as the average spread of a machine should therefore be lower (I'd like to test some more machines, for instance) and I would like to express this in a number. **Questions** Some answer to questions: they aren't identical; if some machine really misbehaves I could see it directly from the machine test (i.e. I would see it from a large X_i), but I wouldn't detect a small misbehavior. Furthermore, the sample amount for each machine could be different (I have more tests on my old machine wrt my new machines).


title: Testing the likelihood that two percentages are equal

text: I'm looking for a method to determine the likelihood that two percentages from independent groups are EQUAL. In this case I'm looking at the likelihood of symptom improvement from patients gathered through two different methods (one traditional double-blind study and one from an online enrollment). So if in the traditional method 60% of 169 patients reported improvement and in the online study 63% of 110 patients reported improvement. Is there a way to tell whether these two percentages are the same? I know I can use a chi- squared test to determine the likelihood that they are different but can I use 1 minus chi-squared as the likelihood that they are the same? Is there a more appropriate test I can do?


title: Can random effects apply only to categorical variables?

text: This questions might sound stupid, but... **is it correct that random effects could apply only to categorical variables** (like individual id, population id, ...), e.g. say $x_i$ is categorical variable: $y_i$ ~ $\beta_{x_i}$ $\beta_{x_i}$ ~ $Norm(\mu, \delta^2)$ but from the principle the **random effect cannot apply to continuous variable** (like height, mass...), say $z_i$: $y_i$ ~ $\alpha + \beta \cdot z_{i}$ because then there is only one coefficient $\beta$ which cannot be constrained? Sounds logical but I wonder why it is never mentioned in statistical literature! Thanks! ~~**EDIT:** But what if I constrain $z_i$ like $z_i$ ~ $Norm(\mu, \delta^2)$? Is it then random effect? But this is different from the constrain I put on $\beta_{x_i}$ - here I constrain the _variable_ whereas in the previous example I constrained the _coefficient_! It starts to look as a big mess to me... Anyway, it doesn't make much sense to put this constraint, because $z_i$ are known values, so maybe this idea is completely odd :-)~~


title: Residuals: marginal or conditional? Raw, Studentized, or Pearson?

text: Residuals are used to check a linear model's assumptions. From SAS Help 'The MIXED Procedure: Residuals and Influence Diagnostics' section we know that there are different types of residuals (Table 58.21 Residual Types Computed by the MIXED Procedure). Which one should we use? What does it mean 'internally' or 'externally' processed? What does 'standardization' mean? What kind of residuals we can calculate with the B formula in my previous question? This, this, this questions, answers, and comments were very helpful but I still don't understand well.


title: OR vs RR when the outcome is common

text: I have a dichotomous outcome (life satisfaction). The proportion of people in my sample who are satisfied is fairly high (about 85%). I know that from my Chi-squared analysis that both income and age are significantly associated with being satisfied (both are categorical variables). I want to run logistic regression with life satisfaction as my outcome and adjust for both income and age to see if they both remain associated with life satisfaction (after adjusting for each other). I'm not interested in the OR estimations per se, only whether they remain associated. I also want to run contrasts to determine differences between groups, post hoc. My main question is: because the outcome is common in the population, should I worry too much about logistic regression inflating the OR compared to the RR? Or should I really be using Poisson regression? My second question is: when I am running contrasts after running logistic regression, should these be on the beta coefficients (the ORs themselves) or should this be between the predicted margins, and will the contrasts be incorrect if the outcome is common? Maybe this doesn't even matter as I don't have a time component...I think I am worrying for nothing


title: Confidence level for a multivariate linear regression

text: I am trying to use the Perl Statistics::Regression package to calculate a multivariate regression. The results I am getting per each of my variables are: `Theta` (the coefficient), `StdErr` (I am not really sure about the meaning of this one) and `T-stat` (which is the coefficient / `StdErr`) On top of that I get `R^2`. Additional parameters the package calculates are: `rsq`, `adjrsq`, `sigmasq`, `set`, and `ybar`. I have two questions:   1. What would signify a good result: high `R^2` value or high `T-stat` values?   2. How do I determine the level of confidence of the result? * * * edit: Thank you very much, chl, for the link below. I think it is very helpful and certainly made several things clearer. I know that my lack of background in statistics is far from ideal but I was hoping I'd still be able to achieve something. Now, if I understand correctly, in order to appreciate the level of confidence it is suggested to examine the F-value (supposed to give an estimation of the model as a whole) and the P-values (supposed to give an estimation of each of the variables). My questions are:   1. Is this correct?   2. What is the exact formula to calculate the P-values? I understand they are dependent on the degrees of freedom but what would be the exact way to calculate it?   3. F = model mean square / error mean square. Is that correct? If so, is the model mean square calculated by dividing the model sum of squares by the model degrees of freedom?   4. I see that the larger the F the more confidence I should have in rejecting the null hypothesis. Any rules of thumb of how large is large? Thanks once again and apologies for the very basic non-intelligent questions...


title: Distribution of distances to an observation from the normal and the median of these distances

text: Given $x\sim \mathcal{N}(\mu=0, \sigma^2=1)$, the squared distances of the $x$ values to $\mu$ are distributed $\chi^2_1$. I am interested in the distribution of the squared distances to an arbitrary $x$. For instance, what is the distribution of the squared distances of the observations to $x=2$? In particular, I would like to know for a given value of $x$ what the median squared distance of the observations to that $x$ is (so, the 50% quantile of that distribution) .


title: Finding a minimum variance unbiased (linear) estimator

text: Here is a basic question that perhaps has a simple answer, but one that I was not able to find by quickly scanning the literature. Suppose that I have a collection of $n$ unopened boxes. Each box $i$ contains either a red ball or a blue ball: write $c_i = 1$ if box $i$ contains a red ball, and $c_i = 0$ otherwise. Moreover, each box $i$ has a probability $p_i > 0$ written on it. If I attempt to open box $i$, then with probability $p_i$ I will succeed and observe whether $c_i = 1$ or $c_i = 0$. With probability $1-p_i$ I will fail, and will not learn $c_i$. I can only try once -- if the box does not open, I never learn $c_i$. The parameters for each box $(c_i, p_i)$ are drawn i.i.d from an uknown distribution $\mathcal{D}$. I wish to produce an unbiased estimate of the expected fraction of red balls:$$C \equiv E_{(c_i, p_i) \sim \mathcal{D}}[ c_i]$$ Note that while the marginal distribution on the $c_i$ is Bernoulli, I make no assumptions about how the $p_i$ are distributed -- they could be correlated with the $c_i$ in some arbitrary way. An obvious first attempt at finding an unbiased estimator is to attempt to open each box, and set $d_i = \frac{c_i}{p_i}$ if the box opens, and $d_i = 0$ otherwise. Then take as my estimate: $$S \equiv \frac{1}{n}\sum_{i=1}^nd_i$$ Since $E[d_i] = C$ for each $i$, this is an unbiased estimator. Can I argue that it is the minimum variance unbiased estimator? If I was literally sampling the values $d_i$, then I could at least say that $S$ was the minimum variance _linear_ unbiased estimator, since each $d_i$ is sampled i.i.d and I am estimating the sample mean. But can I say anything about $S$ given that I am actually sampling $(c_i, p_i)$? Is there a literature on this sort of problem?


title: Why does correlation matrix need to be positive semi-definite and what does it mean to be or not to be positive semi-definite?

text: I have been researching the meaning of positive semi-definite property of correlation or covariance matrix. Any information on   * Definition of positive semi-definiteness    * Its Important properties, practical implications   * The consequence of having negative determinant, impact on multivariate analysis or simulation results etc would be greatly appreciated!