title: How to plug viral coefficient into an exponential growth function?

text: Is there an exponential growth function that has the viral coefficient(k-factor) as a variable? For example I have 10 initial customers Each invites 10 new customers 20% convert into a new customer That gives me a viral coefficient of 2 (10 invites * 20% conversion rate) My 10 customers will send 10 invites each for 100 invites. 20% of those 100 will convert giving me 20 new customers. That leaves me with 30 customers. The exponential growth function is P(t) = P(0)e^kt p(t) is the population at time t p(0) is the initial population k is the growth rate t is the time If I plug in the above numbers I get 30 = 10e^k2 If I solve for k I get = .549 I would have assumed that the growth rate and the viral coefficient would have been the same, but I'm clearly missing something. How does the viral coefficient factor into this equation.


title: Create a custom exponential function

text: I want to create a custom inverse exponential pdf with domain [0.00002, 0.0001] and range [0,1], where x = 0.00002, y = 0 and x = 0.0001, y = 1 The pdf should have the shape like that for a function like 1 - exp(-x). How do I form this function?


title: Back into a viral coefficient from google trends data

text: I want to back into a viral coefficient from Google Trends data. If I look at Foursquare's data, http://www.google.com/trends/explore?q=foursquare#q=foursquare&cmpt=q we see that they went through what looks like an exponential growth curve from July 09 to August 10. I can get the week by week trend score by clicking the gear and then download as CSV. If we make the assumption that all the growth was viral, how would I back into the k-factor using that data?


title: Regression testing after dimension reduction

text: I have a 12 item Likert scale for my predictor variable [IV], and a 9 item Likert scale for my dependent variable [DV]. I used SPSS and did factor analysis on both scales, and found that the IV had 3 components, while the DV had 2. I would like to report the relationship between the IV and the DV using regression analysis. How to accomplish this based of the finds of the factor analysis? Thanks in advance.


title: Need to build a model to see the submersion time

text: I need some help since it is not long time I am using R. I have one year, every day, twice a day (every 6 hours) on tides. I would like to model this in order to check the submersion time of some areas in my lagoon. Could you please tell me how can I do?thank you


title: Creating a model to interpret numerical scores

text: Good morning/afternoon everyone, first of all thanks to all of you for the valuable insights provided. I will be oulining here my current challenge, trying to provide as much detail as possible. **Context:** I have a dataset with scores associated to a list of clients of my company. These scores have been obtained through a survey designed to assess the riskiness of doing buisness with such clients, whereby every answer has a score associated with it. The higher the sum of the scores obtained by each client in each question, the higher the risk profile. Such scores range from 0 (the minimum) to several thousand. **Data quality:** The data I have is not normally distributed and does not fit a normal distribution even if standardized using Zscore and the NORM.DIST function in Excel. Smarter people than me have found distributions that could fit the data (ie: Gamma, Student T, etc.) but these distributions do not transfer over to the dataset of a different time period. That is, the distribution of the data changes every year a survey is taken and I would like to avoid trying to figure out what distribution fits the data everytime a new survey is done. **What I would like to do:** I am interested in creating a way to compare the risk scores associated with each client to come up with a way of deciding what qualified to be high risk and what is low risk. Ideally, I would love to have 5 categorizations: Very High, High, Medium, Low, Very Low risk. I need to devise a way to say that a certain score is high as opposed to low or medium. **What I ruled out:** I ruled out using the arythmectic mean and standard deviation. The reason is that the majority of clients has relatively low scores (0-400) and few have scores in the thousands. I fear that the mean would understate/overstate the mean risk score, being heavily affected by outliers. I need a more robust measure, so I was thinking about using the median instead. **My question:** How do you suggest I go about defining the buckets that I can associate with each of my five desired risk categories (Very High, High, Medium, Low, Very Low risk)? I know I could do it with percentiles relatively simply but that solution is not that appealing to me: it would make somewhat arbitrary to say, for instance, that the top 20 percentiles are very high risk, because it would most likely have a huge gap in between the min and max values of such bucket, whereas this would not be the case with the bottom 20 percentile, since lower score values are much more common. **A note:** I would really like to avoid any solutions that involve attributing a distribution to the data, since that will change and aside from the Gaussian one, most have hard to define parameters, at least in the absense of a statistical software more evolved than Excel. Thank you all very much. Bernardo


title: Order statistics without independence assumption

text: I want to derive an expression for the cdf of the min of a set of $k$ random variables $X_i$ (with same cdf $f_i(x)$) which are identically distributed but not necessarily independent. So far, I got that the cdf of the min is $$ 1−P( \text{all } X_i>x)$$ I know that in the case where they are independent, we can express this as $1−(1−f_i(x))^k$. I was just wondering if we can go further without the assumption that they are ? Thanks !


title: Summarising simulations on a transformed parameter

text: Let $\theta > 0$ be some model parameter for which properties (bias, ...) of an estimate are studied via simulations. For a given data set, an estimate $\hat{\theta}_i$ of $\theta$ can be obtained by maximising a likelihood function. As $\theta$ is positive, however, I perform the maximisation over $\eta = \log(\theta) \in \mathbb{R}$ and I set $\hat{\theta} = \exp(\hat{\eta})$. After running $N$ simulations, how should I compute the mean value,   * $\bar{\hat{\theta}} = \frac{1}{N} \sum_{i=1}^N \hat{\theta}_i$, or   * $\bar{\hat{\theta}} = \exp(\frac{1}{N} \sum_{i=1}^N \hat{\eta}_i)$? * * * **Example** : \-- Minus log likelihood for a sample from the Exp distribution with mean $\theta = \exp(\eta)$:               minusloglik <- function(eta, sample)     {       theta <- exp(eta)       - sum(dexp(x=sample, rate=theta, log=TRUE))     }      \-- True value of $\theta$:               theta <- 5.73      \-- Simulations:               thetaHat <- etaHat <- rep(NA, 1000)     for(i in 1:1000)     {         sample <- rexp(n=100, rate=theta)       etaHat[i] <- nlm(f=minusloglik, p=0, sample=sample)$estimate       thetaHat[i] <- exp(etaHat[i])     }      **Question** :   * Should I summarise the results as `mean(thetaHat)` or as `exp(mean(etaHat))`?   * Is the answer the same if $\theta$ denotes the variance of a normal distribution?


title: If three random variables have the same variance, what will the co-variances look like?

text: I'm curious to know, if three random variables have the same variance, what will the co-variances look like? Can somebody help me to figure it out?


title: What am I really calculating here?

text: This is partly related to the question at Visualising the variance which I haven't accepted an answer yet, partly because I have become unsure if I'm asking the correct question. The case is this, we have a population P that consists of 50 people. These have answered questions on a scale of 1-5, which our software normalizes to a 0-100 scale. We also have a reference population R which consists of several hundred, if not thousands of people. The reference data is "pre-rendered" to have an average score, a standard deviation, and an average standard deviation (ASD). The last value is found by doing the standard deviation per survey, then the average of those standard deviations. (The reference population contains a representative collection of previous surveys) We then generate a value for what we directly translated call "variation", but I have started to believe is more a comparison of spread... This is where my issues lie. The formula related to this is               (((((stdev(P) - ASD) / ASD) * 100) + 100) / 2)      (The last part is to adapt to a 0-100 scale, where 50 means that the spread in your result is the same as the spread in the reference population) What are we really calculating with this formula? (I have taken over this system as a developer, and statistics are not my strong side)


title: Is spurious regression impossible if you include lagged error terms?

text: If you have a random walk A, and a random walk B, and you regress them against each other you run into spurious regression. In our textbook however is written that it is impossible to have spurious regression in dynamic models. So in case you regress A against lagged terms of A plus B and lagged terms of B, you don't have the problem of spurious regression?


title: Adding quadratic term changes the sign of the variable

text: Number of books published in a year (noBook) is my dependent variable and I have independent variables including the age of the author (age). The coefficient of age is positive and it is significant. When I add age^2 to the model, age and age^2 are significant but the sign of age becomes negative. Could you please advise why it is like that? And, which model should I consider?


title: Parameters’ uncertainty for small sample size

text: Suppose we have a small set of numbers (5 to 10 observations), and we’re trying to fit a distribution to this set. Also, we know that all numbers are positive. I tried to fit lognormal, but I’m not sure how good my estimates are since the sample is very small; also I’m not sure how whether it is enough to look at goodness-of-fit test due to the small sample size. Any suggestions on how to tackle this issue (i.e., to be more confident, certain, about my estimates)?


title: clustering evaluation based on a gold standard

text: I have a dataset of items and I want to measure how my clustering method works. I'm using R and simple k-means clustering. I have a lets say gold standard of my clusters and I want to see how good is the result of my clustering algorithm based on the features I used. Is there any straightforward way to do that in R? I was thinking about something like Jaccard similarity per item. I mean the intersection/union of clusters containing each item and somehow mixing them together!


title: Mixing probabilities in mixture models using EM

text: Assume we want to estimate the mixing probabilities ($\pi_{k}$) for each member distribution in the mixture model. We know that $\sum_{m}^{K}\pi_{m}=1$, so we can formulate the optimization problem as a Lagrangian multiplier; the function we want to maximize (w.r.t. $\pi_{k}$) is is given by, $$ Q = \sum_{l}^{n}\sum_{m}^{k} p(m|l)\text{ log}q(m|l)+\delta(\sum_{k}\pi_{k}-1) $$ where: $q(k|l) = \pi_{k}g_{k}(x_{l};\theta_{k})$, i.e. the joint probability of selecting mixture component $k$ and the probability of datapoint $x_{l}$ given component distribution $g_{k}$ and $p(k|l) = \frac{q(k|l)}{\sum_{m}^{k}q(m|l)}$ Setting $\frac{\partial Q}{\partial \pi_{k}}=0$, yields, $$ \frac{\partial Q}{\partial \pi_{k}} = \sum_{l}^{n}\frac{p(k|l)}{\pi_{k}}+\delta = 0 $$ Now, in order to solve for $\pi_{k}$, we must solve for $\delta$; apparently, $\delta=-n$. This results in, $$ \pi_{k} = \sum_{l}^{n} \frac{p(k|l)}{n}$$ However, I'm not sure why this is; more specifically, how can it be shown that, $$ \sum_{l}^{n}\frac{p(k|l)}{\pi_{k}} = n $$


title: Interpretation of regression analysis with integer value as response and predictor

text: I have a data set of partygoers. For each partygoer I know which party she attended. I want to be able to estimate how many partygoers a pair of parties ($i, j$) _should_ have in common _given_ the number of participants of each of the two party. The idea is to spot within the dataset pair of parties which have higher or lower shared partygoers than the average controlling for the size of both parties. Let's be clear. I _do_ know who went to which party. I _don't_ want to predict something that is not in the dataset. What I _do want_ to know is whether the number of shared participant by each pair of parties is far from the "average" of the dataset. But somehow I need to control for the size of the two parties because of course the higher the size the higher the chance to see people being in both parties. I imagined a linear regression model like PartygoersInCommon $\approx \beta_0 + \beta_1 ($Partygoers$_i$ Partygoers$_j)$ I have two questions concerning the interpretation of the results: 1) The fact that the response PartygoersInCommon is not a continuous variable but an integer does complicate the interpretation of the model? 2) Is it correct to multiply Partygoers$_i$ and Partygoers$_j$ or would be better to treat them as two separate independent variables like in PartygoersInCommon $\approx \beta_0 + \beta_1$ Partygoers$_i + \beta_2$ Partygoers$_j$?


title: Order of preprocessing steps in a binary classification problem

text: I have these stages (ordered) for `preprocessing` in my binary classification problem.   1. Dividing data based on criteria (class1 and class2 databases)   2. Outlier detection `(removing X - mean(X) >= 2*std(X) based on every feature)`   3. data balancing ( in every year - I have financial data and I'm balancing data based on some features of samples - database consisted of 5 years and I'm doing this step separately in every year - now i have equal samples in every class)   4. Data normalization   5. Dimension reduction (`PCA` etc.) Is this ordering true? You know that `normalization`, `dimension reduction` and `outlier detection` generate some parameters that we must apply these parameters to new database that will import to our created database in future( out sample in real using of the model raw data without any balancing) . I know that `outlier detection` should be before `normalization`. Is this two `preprocessing` steps (`outlier detection` and `normalization`) should be before `balancing` of after `balancing` of data?


title: Computing repeatability of effects from an lmer model

text: I just came across this paper, which describes how to compute the repeatability (a.k.a. reliability, a.k.a. intraclass correlation) of a measurement via mixed effects modelling. The R code would be:               #fit the model     fit = lmer(dv~(1|unit),data=my_data)          #obtain the variance estimates     vc = VarCorr(fit)     residual_var = attr(vc,'sc')^2     intercept_var = attr(vc$id,'stddev')[1]^2          #compute the unadjusted repeatability     R = intercept_var/(intercept_var+residual_var)          #compute n0, the repeatability adjustment     n = as.data.frame(table(my_data$unit))         k = nrow(n)         N = sum(n$Freq)     n0 = (N-(sum(n$Freq^2)/N))/(k-1)          #compute the adjusted repeatability     Rn = R/(R+(1-R)/n0)      I believe that this approach can also be used to compute the reliability of effects (i.e. sum contrast effect of a variable with 2 levels), as in:               #make sure the effect variable has sum contrasts     contrasts(my_data$iv) = contr.sum          #fit the model     fit = lmer(dv~(iv|unit)+iv,data=my_data)          #obtain the variance estimates     vc = VarCorr(fit)     residual_var = attr(vc,'sc')^2     effect_var = attr(vc$id,'stddev')[2]^2          #compute the unadjusted repeatability     R = effect_var/(effect_var+residual_var)          #compute n0, the repeatability adjustment     n = as.data.frame(table(my_data$unit,my_data$iv))     k = nrow(n)     N = sum(n$Freq)         n0 = (N-(sum(n$Freq^2)/N))/(k-1)          #compute the adjusted repeatability     Rn = R/(R+(1-R)/n0)      Three questions:   1. Do the above computations for obtaining the point estimate of the repeatability of an effect make sense?   2. When I have multiple variables whose repeatability I want to estimate, adding them all to the same fit (e.g. `lmer(dv~(iv1+iv2|unit)+iv1+iv2`) seems to yield higher repeatability estimates than creating a separate model for each effect. This makes sense computationally to me, as inclusion of multiple effects will tend to decrease the residual variance, but I'm not positive that the resulting repeatability estimates are valid. Are they?   3. The above cited paper suggests that likelihood profiling might help me obtain confidence intervals for the repeatability estimates, but so far as I can tell, `confint(profile(fit))` only provides intervals for the intercept and effect variances, whereas I would additionally need the interval for the residual variance to compute the interval for the repeatability, no?


title: What university level statistics courses are considered advanced/hard?

text: I've recently been looking for top-of-the-line statisticians in a recruiting process for our company. Myself, I'm a Physics Engineering major. I gather that great mathematical statisticians have studied a bit different courses, and much more in depth. When evaluating a candidate, are courses a good indicators of this person being excellent? Preferably we're talking graduate or post-graduate level. * * * We're looking to fill roles of data miners, statistical modeling and data visualization. Thanks Chris, for the suggestion to clarify.


title: Exponential weighted moving skewness/kurtosis

text: There are well-known on-line formulas for computing exponentially weighted moving averages and standard deviations of a process $(x_n)_{n=0,1,2,\dots}$. For the mean, $\mu_n = (1-\alpha) \mu_{n-1} + \alpha x_n$ and for the variance $\sigma_n^2 = (1-\alpha) \sigma_{n-1}^2 + \alpha(x_n - \mu_{n-1})(x_n - \mu_n)$ from which you can compute the standard deviation. Are there similar formulas for on-line computation of exponential weighted third- and fourth-central moments? My intuition is that they should take the form $M_{3,n} = (1-\alpha) M_{3,n-1} + \alpha f(x_n,\mu_n,\mu_{n-1},S_n,S_{n-1})$ and $M_{4,n} = (1-\alpha) M_{4,n-1} + \alpha f(x_n,\mu_n,\mu_{n-1},S_n,S_{n-1},M_{3,n},M_{3,n-1})$ from which you could compute the skewness $\gamma_n = M_{3,n} / \sigma_n^3$ and the kurtosis $k_n = M_{4,n}/\sigma_n^4$ but I've not been able to find simple, closed-form expression for the functions $f$ and $g$. * * * **Edit:** Some more information. The updating formula for moving variance is a special case of the formula for the exponential weighted moving covariance, which can be computed via $C_n(x,y) = (1-\alpha) C_{n-1}(x,y) + \alpha (x_n - \bar{x}_n) (y_n - \bar{y}_{n-1})$ where $\bar{x}_n$ and $\bar{y}_n$ are the exponential moving means of $x$ and $y$. The asymmetry between $x$ and $y$ is illusory, and disappears when you notice that $y-\bar{y}_n = (1-\alpha) (y-\bar{y}_{n-1})$. Formulas like this can be computed by writing the central moment as an expectation $E_n(\cdot)$, where weights in the expectation are understood to be exponential, and using the fact that for any function $f(x)$ we have $E_n(f(x)) = \alpha f(x_n) + (1-\alpha) E_{n-1}(f(x))$ It's easy to derive the updating formulas for the mean and variance using this relation, but it's proving to be more tricky for the third and fourth central moments.


title: Where to find mathematical modeling help on low-budget project?

text: (Not sure how to formulate this question and also didn't find any suitable tags.) Are there any online sites where one can find professionals who do mathematical modeling, analyzing data and so on? More general, where to turn if you need help modeling some data on a (low budget) project?


title: Censoring/Truncation in JAGS

text: I have a question on how to fit a censoring problem in JAGS. I observe a bivariate mixture normal where the X values have measurement error. I would like to model the true underlying 'means' of the observed censored values. \begin{align*} \lceil x_{true}+\epsilon \rceil = x_{observed} \ \epsilon \sim N(0,sd=.5) \end{align*} Here is what I have now:                for (i in 1:n){        x[i,1:2]~dmnorm(mu[z[i],1:2], tau[z[i],1:2,1:2])        z[i]~dcat(prob[ ])      }      Y also has measurement error. What I want to do is something like this:                for (i in 1:n){        x_obs[i] ~ dnorm(x_true[i],prec_x)I(x_true[i],)        y_obs[i] ~ dnorm(y_true[i],prec_y)        c(x_true[i]:y_true[i])~dmnorm(mu[ z [ i ],1:2], tau[z[i],1:2,1:2])        z[i]~dcat(prob[ ])      }           #priors for measurement error      e_x~dunif(.1,.9)      prec_x<-1/pow(e_x,2)      e_y~dunif(2,4)      prec_y<-1/pow(e_y,2)      Obviously the c command is not valid in JAGS. Thanks in advance.


title: Soft Margin Loss and Conditional Probabilities

text: The following question is taken almost verbatim from "Learning with Kernels", by Scholkopf and Smola. Consider the following soft margin loss function: $max(0,1−yf(x))$. I have the problem of needing to compute the conditional probability p(y|x) corresponding to this function and am having trouble making the connection between conditional probability and this function. I have come across a few papers (e.g. http://www.unc.edu/~yfliu/papers/lum.pdf) that say that "soft classifiers explicitly calculate the class conditional probabilities", but I do not understand how as the output of a classifier is not a probability. Can someone please explain what I am missing?


title: Bias in EM estimation for a mixture of normal distributions

text: Are the parameter estimates for a mixture of two normal distribution using EM algorithm biased or unbiased? More specifically, if I use the EM-algorithm to obtain ML estimates of $μ_1$, $μ_2$, $σ^2_1$, $σ^2_2$ and $p$ in a two component normal mixture model, are the ML estimates of those five parameters unbiased or are some or all of them biased? If any are biased, how do I compute the bias?


title: GLS, heteroskedasticity and Ridge Regression/Lasso

text: I am hoping to use a regularised regression technique, using cross validation, to fit a linear model to a set of predictors which have some highly correlated variables. However, I also know (highly suspect) that there is a significant amount of heteroskedasticity in the errors. I saw this question: Ridge and LASSO given a covariance structure?, which suggests that if I know the variance structure it is easy to incorporate it into these procedures. However I do not, and instead was hoping to estimate it using a GLS or Feasible-GLS framework. My initial thought is to run my CV-regression process with no variance model, estimate one from the residuals, scale and re-run the regression. However this feels very ad-hoc and I'd like to know if there is an obvious problem with this approach?


title: Question about Neural Network Training

text: **What is the benefit of optimizing our neural network error function using back prop rather than just using gradient or stochastic gradient descent directly on the error function?** How come we don't just take the partial derivative of the error function with respect to every weight and do gradient or stochastic gradient descent based on that? Instead of doing back propagation to come up with these artificial error values for hidden nodes and then for each hidden node multiplying its error times the various partial derivatives of weights associated to its incoming edges, as if that hidden node was actually an output node.


title: Weighting, standardize, transforming data before for Multivariate ANCOVA

text: My question is in regards to data manipulation prior to Multivariate ANCOVA. My research focuses on crustacean eggs. The model involves a covariate and 2 x independent variables and 4 dependent variables. The dependent variables are to do with egg traits. Two of the continuous dependent variables (DV) are known to increase with egg stage (1-4)(not the covariate and not one of the DVs). DV 1 is egg size, the egg samples have approx 25 eggs measured for size and the sample is assessed for egg stage. lets say 15 of the 25 eggs a stage 2 and 10 are stage 3. so a proportion of the sample is stage 2 and another proportion stage 3, I have mean size for each stage taken from the entire data set. I want to reduce the effect egg stage on size or standardize the diameter for egg size some how before putting the data in the MANCOVA. Your help is appreciated. Thanks


title: Some doubt in reading Machine Learning A Probabilistic Perspective ( chapter 3.2 )

text: When I am reading Murphy's Machine Learning A Probabilistic Perspective. In chapter 3.2. I have some doubt. I think the author want to express is two things. First, we can use Bayes formula to computer posterior, combine likelihood and prior. Second, if data is big enough, we can compute it directly. $1$ or $0$ depend on the hypothesis which maximum itself($p(h|\mathcal{D})$) is exist or not. I can't see these conclusion clearly. And I feel indigestion. the below is some detail I can't go through. * Likelihood. He give an formula:$$p(\mathcal{D}|h)=\left[\frac{1}{\text{size}(h)}\right]^N=\left[\frac{1}{|h|}\right]^N$$ And the author say this is the probability of independently sampling $N$ items(with replacement) from $h$. Yes, if the probability mass function is uniform. this is equal to $\left[\frac{1}{|h|}\right]^N$. but I think the the notation $p(\mathcal{D}|h)$ is not reasonable. It represents the probability of event $\mathcal{D}$ when $h$ is happen. under the assumption **pmf** is uniform. it must equal to $$\frac{\text{size}(\mathcal{D})}{\text{size}(h)}$$. * Posterior , the formula is $$p(h|\mathcal{D})=\frac{p(\mathcal{D}|h)p(h)}{p(\mathcal{D})}$$ It is exactly Bayes formula. Then ,the author say if we have enough data, the prior is not important. the likelihood almost determine everything.then he give a formula: $$p(h|\mathcal{D})\rightarrow \delta_{\hat{h}^{MAP}}(h)$$ where $\hat{h}^{MAP} = \max_{h}p(h|\mathcal{D})$ and $\delta$ is Dirac measure. so does it mean if exist hypothesis $\hat{h}$ in $\\{h\\}$ maximum $p(h|\mathcal{D})$. then $p(h|\mathcal{D}) \rightarrow 1$? To find $\hat{h}$ maximum $p(h|\mathcal{D})$ is equal to maximum $\log p(\mathcal{D}|h) + \log h$. Then the author say if we have enough data, only to find $\hat{h}$ maximum $\log p(\mathcal{D}|h)$ is enough. Since the likelihood term depends exponentially on $N$, and the prior stays constant. But I don't think so. Intuitively, If we have a lot of data. the prior is not important. But for this formula, both $p(\mathcal{D}|h)$ and $p(h)$ are probability. they $\in[0,1]$. Although $\mathcal{D}$ is big. for different hypothesis $h$, $p(\mathcal{D}|h)$ is different. but it $\in [0,1]$ the same as $p(h)$. So I think we can't ignore $p(h)$.


title: What is the difference between a multi-label and a multi-class classification?

text: What is the difference between multi-label classification and multiclass classfication. Speficially, what is the difference between a label and a class? Please provide a clear example. _"Multiclass classification should not be confused with multi-label classification, where multiple labels are to be predicted for each instance." -wikipedia_ ... was not very helpful.


title: Testing for significance between means, having one normal distributed sample and one non normal distributed

text: I have following problem: Within an independent groups 1-factor design I have two independent groups, with a sample size of 20 each. The data of the treatment group is not normally distributed, whereas the data for the control group is (checked with Shapiro- Wilk Normality Test). Now I want to check if the differences of the means of both groups are significant. What is the appropriate test for this? I think it should be the Wilcoxon Rank Sum and Signed Rank Test, but I am not sure... Could please anybody help me?


title: On what tasks does neuroevolution outperform basic application of neural networks or genetic algorithms?

text: There has been a recent interest in combining genetic algorithms and neural networks into a general neuroevolution framework. The basic idea, is that your genetic algorithm is evolving the parameters of many neural-network which are then used to solve your task at hand. A sort of genetic programming but instead of evolving a snippet of code to do some task, you are evolving a neural network. When should I use this combined approach instead of using neural-networks or genetic algorithms by themselves? For what types of problems has the combined approach produced better results that the individual approaches? For what types of problems is the combined approach the 'best' approach?


title: Interpreting p-values associated with correlation measurements

text: The "Introductory Statistics with R" book contains a section that deals with correlations (section 6.4 in the second edition). The book shows Pearson, Spearman and Kendall correlation coefficients computed on the `blood.glucose` and `short.velocity` columns of the thuesen data set. The p-values associated with these coefficients are 0.048, 0.139 and 0.119, correspondingly. The book then says the following: > Notice that neither of the two nonparametric correlations is significant at > the 5% level, which the Pearson correlation is, albeit only borderline > significant. I have several problems with this paragraph. First of all, my naive guess would be that since the non-parametric coefficients do not imply linearity, they will tend to be "significant" more frequently than Pearson's r. Am I right? Secondly, and more importantly, is such a comparison between p-values of different tests applied on the same data legit? (I'm talking about real-life comparisons and not about trivial examples in a text book) If it is, how one need to interpret the notion that linear correlation is "significant", while rank or concordance correlation isn't?


title: Testing significant difference between two determination coefficients $R^{2}$ for two PLS-based structural equation models

text: Two structural equation models were tested (one was based on a sample with 199 individuals and the second one on a sample with 93 individuals). The aim was to test whether the results of the first model could be replicated using a new dataset. **Is there a way to calculate whether the $R^2$ (of the dependent variable) differ significantly between both models?** Is it possible to use Fisher Transformation for this purpose? Using the equations behind this Tool I can use R to calculate z and a corresponding p-value: $z_{A} = \frac{1}{2} \ln\left({1+\sqrt{R^2_A} \over 1-\sqrt{R^2_A}}\right) = \operatorname{artanh}(\sqrt{R^2_A})$ $se(z) = \sqrt{\frac{1}{n_a - 3} + \frac{1}{n_b - 3}}$ $z = \frac{z_A - z_B}{se(z)} \tilde{} t(n_A + n_B - 2)$               fisher_z <- function(rA, rB, nA, nB){       z_diff <- (atanh(sqrt(rA)) - atanh(sqrt(rB))) / sqrt((1/(nA-3)) + (1/(nB-3)));       pvalue <- pt(z_diff, nA + nB - 2);       return(list("z_diff"=z_diff, "p"=pvalue))     }


title: Including a baseline measure of the DV as a covariate in an ANCOVA

text: I have 4 groups based on BPD symptoms at baseline and follow-up. (BPD T1 and T2; BPD T1 not T2; no BPD T1 but BPD T2; Never BPD). I am examining whether psychosocial functioning outcomes (SOFAS score) differ amongst the groups. As I wish to control for Axis I disorders (mood, anxiety, and substance use) and treatment (service 1 or 2) I am using an ANCOVA. I would also like to control for the baseline measure of the outcome variable (SOFAS at T1). The Axis I disorders and the treatment are all dichotomous variables; SOFAS scores (both T1 and T2) are continuous. When using SPSS is it correct to put the dichotomous variables in the "fixed factors"box? ANCOVA Fixed factors = Groups; mood; anxiety; substance use; service Covariate = SOFAS T1 DV = SOFAS T2. Does this seem reasonable? From my readings it seems OK but my supervisor is unsure and wants me to get further advice. So I would appreciate any comments. Thank you


title: What does Jaynes' continous pdf notation "g(x)dx" actually mean?

text: Something has been bugging me about E.T. Jaynes' treatment of continuous parameters. In his book Probability Theory: The Logic of Science, uses notation that I am unfamiliar with when getting probabilities for continuous parameters from their densities. For example, in chapter 4, page 420, equation 4-49, Jaynes "discretizes" the continuous parameter $f \in [0,1]$ by taking the interval $(f, f+df)$ and putting a probability on it via $g$, the density of $f$: > .. so we will go back to the original probability from (4-3): > > $P(A|DX)=\frac{P(D|AX)P(A|X)}{P(D|X)}$ > > Letting A now stand for the proposition "The fraction of bad ones is in the > range $(f, f+df)$", there is a prior pdf > > $P(A|X)=g(f|X)df$, > > which gives the probability that the fraction of bad ones is in the range > $df$, Why isn't the above written as $P(A|X)=\int_f^{f+dx} g(t|X)dt$ (with dummy variable $t$)? $g(f|X)df$ for a finite $df$ is the left-rectangle approximation (a la Riemann sums), not $P(A|X)$. If you were to partition the $[0,1]$ into many separate proposition and assigned probability as above, it may not sum to 1. Or does $df$ denote some sort of infinitesimal notation I'm not familiar with? I've never seem that notation used 'alone'-- afaik it's only meaningful when used in context of differentiation/integration, like $\frac{d}{dx}$ or $\int foo()dx$. If it is some infinitesmal notation, this seems to run contrary to Jaynes' philosophy of sticking to finite number of propositions until the end, where the limit may be taken explicitely. As another example, in chapter 15, page 1514, equation (15-38), he writes the bivariate normal probability with correlation $\rho$ as: > $p(dx\ dy|I)=\frac{\sqrt{1-\rho^2}}{2\pi}exp[-\frac{1}{2}(x^2 + y^2 -2 \rho > x y)]dxdy$ Which is a little different again. Presumably $p(dx\ dy|I)$ means "Probability that the true value of $(x,y)$ falls in $(x+dx,y+dy)$", but here again he writes the Riemann-sum-like approximation instead of the integral of the density on $[(x,x+dx) \times (y,y+dy)]$. What am I missing?


title: Measuring the smoothness of time series

text: I work on a method that gives a (noisy) estimation of brain volume over time in Alzheimer's patients. As we know that the evolution is smooth and even mostly linear if looked at over a time frame of a few years, one way of evaluating the algorithm is to look at the smoothness or linearity of the estimated time series (brain volume over time). This is why I'm looking for a good metric of smoothness or linearity of time series. What I've been using for now is the average R² from patient-wise linear regressions. This does reflect linearity, but it depends heavily on the sampling (it will give better results with fewer data points). Could anyone suggest a better metric? **Some details since apparently my question wasn't clear**. I have an algorithm that gives a few (from one to five or six) discrete $V(t_n)$ brain volume measurements for many patients. I want to measure to which degree $V$ is close to something linear, eg. $V(t_n) = a t_n + b$ accross all patients (all patients will have different coefficients, but I want an overall measure of linearity). I should also mention that the $t_n$ values are not regularly spaced. **One more comment** Obviously patients with one or two data points are not informative. However I want to take advantage of all the other patients. Even a patient with three data points is informative: three aligned points are much better than three completely jagged points.


title: Fitting a beta-binomial model in the case of overdispersion in R

text: I'm estimating some count data. I have counts for say $m=100$ individuals. Unfortunately when using the Poisson regression overdispersion occurs. So I was thinking to fit a negbin model. But this is not appropriate in my case. So I assume that I can not fit a Poisson regression, because the way the Poisson distribution arises is not appropriate in my case ($n$ is not growing to infinity and $p$ is not converging to zero). So I found the beta-bin model. But quite honestly I'm absolutely not familiar in estimating beta-binomial models using R? First of all: Does it make sense to fit a beta-bin model when anyone wants to estimate counts? Btw: If it makes sense, does anybody know a good book where the application is discribed?


title: Silhouette coefficients after deleting some data and re-clustering

text: I was working on a dataset today on which I used K-means clustering algorithm and then calculated the Silhouette coefficients for each point. I then removed 5% of the data with worst silhouette coefficients and re-clustered again. The average silhouette coefficient for each cluster got worse after re-clustering and removal of those outliers (I expected opposite). Am I getting the right answer or have I made a blatant mistake? The results I got don't seem intuitive to me.


title: self-selection bias due to nonresponse?

text: I have administrative data from the whole population of new doctorates, in a given year, from my region. We have also survey data from a sample of this same population (where the whole population was contacted to participate in the survey, but only around 65% of doctorates participated: this is the sample I should work on). There will be self-selection bias if the non-response is not random. I want to investigate possible wage gaps between different groups in my region (the dependent variable is the log of the ratio between two average wages). Any inputs on how to tackle this issue? Literature on this issue?


title: Interpreting a plot of n confidence intervals for the mean and variance

text: I wanted to know how I would go about interpreting a plot such as this: ![Confidence Interval Plot](http://i.stack.imgur.com/cCrqr.png) where the vertical line segments represent 50 realizations of a confidence interval for $\mu$ How would I explain the location of $\mu$ with regards to the line segments? Assume the lines were more to the left or right of $\mu$, what would that signify?


title: Effect of smaller sample size on confidence intervals

text: Assume we were to collect 20 data samples for the passage time of light and in those 20 samples, there were 3 negative values. If we were to discard these negative values, what effect would that have when constructing confidence intervals for the mean and variation? I know that a smaller sample size results in a larger confidence interval but what if the data is flat out incorrect (i.e. negative times such as in this example)?


title: Simulating data for mixed models

text: Simulating data for linear regression is fairly easy. You generate a vector of X values of some range, you decide on some beta value and you generate random errors, and build a linear equation. Now I want to do the same for a mixed model, with a random effect, and repeated measures. Why I want to do this ? For sample size estimation, I want to check how large should N be in my coming study. I do not know how to create the simulation. I don't mind using R or SAS, whatever you guys can help me with, fine with me. Some words on my design. I have a continuous response variable Y. I have a fixed factor with two levels, I will call it X1 (treatment vs. control). These treatments are some operating techniques. There are two kind of operations these treatments can be applied on, A and B. This is a random factor, because there could be more, but only 2 will be examined. I will call the random factor Z1 (A vs. B). Each subject will have 2 operating procedures (repeated measures) of either type A or B (both repeated measures from the SAME type - both using SAME treatment, i.e., the new treatment or control). I assume a correlation of 0.8 between repeated measures of the same subject, for all subjects. Let's also say I know or predict the means and standard deviations of the groups of X1, and the variance of Z1. How do I simulate or generate data to fit this description so I can run a mixed model on it using various sample sizes and calculate the power ? Thank you ! P.S My main problem is generating the correlations and to take into account the random effect.


title: Follow up test for 3 way interaction

text: This question seems to have been asked before in different ways, but the answers are all not very helpful. In general, suppose three are three factors ($A$, $B$, and $C$), each with two levels. The three way interaction is significant. Then there are $6$ possible two way interactions (perhaps combinations of those interactions?) that could explain the three-way interaction.               A1: BxC,     A2: BxC,     B1: AxC,     B2: AxC,     etc...      What is the proper way to explore these simple interactions? Just run the interactions, run the interactions with some correction for multiple comparison, or calculate the $F$ using some other error term in the denominator? In general, one would think that if $B\times C$ were significant at $A_1$ but not $A_2$, then this would be the best explanation of the three way interaction. . . right?


title: Assessing whether 2 models differ "substantially" (though equal RMSE), to improve prediction

text: In preparation of stacking several models, I would like to gather some heterogeneous models. I am not sure, but expect that the stacking will improve when I have very different models, though all have a (relatively) small RMSE. Imagine that 2 models produce the fitted values on the training set, as the truth is known:               Model1  Model2   Reality       10      5         5        6      3         3        6      3         6       10      5        10      Below the RMSE per model are given. However, I would like some measure of how much the models themselves are similar (here I thought it is possible to do a RMSE between them)               Model1 RMSE: 2.915  (RMSE(model1, reality))     Model2 RMSE: 2.915  (RMSE(model2, reality))     Model1 vs Model2: 4.123 (RMSE(model1,model2))      I would like to find some way to express the difference in terms of how relevant it is to the current RMSE that I try to minimize. The point being that when models are too similar then there is no point in including them all in the stacking, while if there is variation, it might be worth to stack. So, in what way is it possible to assess whether models differ "enough"; is there some formula?


title: Differentiating correlation and causation using conditional probablity

text: I'm trying to understand the difference between causation and correlation using conditional probabilities. From what I understand, one may quantify causation by $P(E_1|E_2) / P(E_1)$. For example, if this ratio is $> 1$, then $E_2$ increase the likehood of $E_1$. Now I try to come up with something similar explaining correlation and I thought of : $$P\big(\left(E_1 \cap E_2\right) \cup \left(E_1^c \cap E_2^c\right)\big)$$ where $E_1^c$ is the complement of $E_1$. But I can't think of correlation in terms of conditional probability. Does correlation alone give information about conditional probability ?   Thanks for helping.


title: How can I align/synchronize two signals?

text: I'm doing some research into respiratory physiology but have come stuck at the analysis stage (should've paid more attention to my stats lectures). I've collected two simultaneous measures of ventilation: the volume of air passing through the mouth and the change in resistance of a chest band to measure chest expansion. I'd like compare the signals and ultimately hope to derive volume from the chest expansion signal. But first I have to align/synchronise my data. As recording doesn't start at precisely the same time and chest expansion is captured for longer periods I need to find the data that corresponds to my volume data within the chest expansion data set. While I could eyeball it, I'll be having to repeat this process god-knows how many times and would like to have a measure of how well they are aligned. I'm not quite sure how to go about this without an anchorpoint (eg. same start time) and with the data in different scales and at different resolutions. I've attached an example of the two signals (https://docs.google.com/spreadsheet/ccc?key=0As4oZTKp4RZ3dFRKaktYWEhZLXlFbFVKNmllbGVXNHc), please let me know if there's anything further I could provide. Apologies if this question is a little vague or convoluted, it's been a long, long day and my words are failing me. I'll attempt to clean it up tomorrow morning, but would welcome any suggestions you had in the meantime. Thanks very much!


title: Assign weights to variables in cluster analysis

text: I want to assign different weights to the variables in my cluster analysis, but my program (Stata) doesn't seem to have an option for this, so I need to do it manually. Imagine 4 variables A, B, C, D. The weights for those variables should be               w(A)=50%     w(B)=25%     w(C)=10%     w(D)=15%      I am wondering whether one of the following two approaches would actually do the trick:   1. First I standardize all variables (e.g. by their range). Then I multiply each standardized variable with their weight. Then do the cluster analysis.   2. I multiply all variables with their weight and standardize them afterwards. Then do the cluster analysis. Or are both ideas complete nonsense? **[EDIT]** The clustering algorithms (I try 3 different) I wish to use are k-means, weighted-average linkage and average-linkage. I plan to use weighted- average linkage to determine a good number of clusters which I plug into k-means afterwards.


title: Finding optimal cutoff point with two linear regression models

text: I am trying to get an optimal cut-off value dividing group with minimum sums of squares of residuals (=observed y - estimated y) the model is like below. > In group 1 : model y= a1x + b1z + C1v ... > > In group 2 : model y= a2x + b1z + C1v ... I have the data of y, x, z, v... The problem is the group 1 and 2 are not divided yet and the purpose of the analysis is finding optimal cut-off point of x using regression models. I searched again and again, but couldn't find the way to make models varying 'a' and share b1 and c1... and fitting it to data. I asked similar quesion in stackexchange, and somebody advised me the problems of this kind of approach, however, I need this approach, because it's some clinical research want to 'find' optimal (not perfect) cut-off point of x. The article I read described below The authors of the article mentioned that they used R, but I cannot find any reference or examples about this kind of analysis. > To determine the relationship between serum 25(OH)D and iPTH concen- > trations while adjusting for confounders that could affect serum 25(OH)D > concentrations (i.e., age, gender, body weight, calcium intake, physical > activity, and season of year), we considered two linear regression models, > one for subjects below a certain concentration of serum 25(OH)D and the > other for subjects above that concentration. To determine the specific > cutoffs, we fitted the two linear regression models described above and > calculated the sums of squares of residuals (=observed PTH - estimated PTH) > from the two models for each concentration of serum 25(OH)D. The models with > the lowest residual sums of squares were our best models, and the > corresponding concentrations of serum 25(OH)D were defined as the optimal > cutoff values. Somebody said that this question is already answered in "regression model fitting for define cut-off" but, I don't think so... It's not regression discontinued design, because there is no a-priori cut-off. Finding cutoff is the purpose of analysis. Thanks.


title: How to measure the influence of an index on the correlation between 2 other indices

text: I have 3 raster datasets. All 3 are values of continuous indices varying between 1 and -1. These are satellite measurements of the same area, divided in a raster. Each cell has a value for every index. The indices representing measures in a pixel are:   * The `NDVI`, representing vegetation.   * The `NDWI`, representing water.   * The `SWI`, representing temperature. Now I would like to know the influence of `NDWI` (water) on the correlation between `SWI` (temperature) and `NDVI` (vegetation). Temperature has a positive relation with vegetation, water a negative relation with temperature. Water has a positive relation with vegetation. Now I would like to show that water has a negative effect on the correlation between vegetation and temperature. So the temperature will go down and vegetation will go up when introducing a higher water value. So more water will show a lower value of vegetation/temperature correlation. Could I just correlate `NDWI` (water) with the correlation results of `SWI` (temperature) and `NDVI` (vegetation)? I have looked at the Fisher z test, but this is just about comparing two correlations, and so I think this does not cover my question. I hope you could give me some advice.


title: Simulating data to fit a mediation model

text: I am interested in finding a procedure to simulate data that are consistent with a specified mediation model. According to the general linear structural equation model framework for testing mediation models first outlined by Barron and Kenny (1986) and described elsewhere such as Judd, Yzerbyt, & Muller (2013), mediation models for outcome $Y$, mediator $med$, and predictor $X$ and are governed by the following three regression equations:   1. $Y = b_{11} + b_{12}*X + e_1$   2. $med = b_{21} + b_{22}*X + e_2$   3. $Y = b_{31} + b_{32}*X + b_{32} * med + e_3$ The indirect effect or mediation effect of $X$ on $Y$ through $med$ can either be defined as $b_{22}*b_{32}$ or, equivalently, as $b_{12}-b_{32}$. Under the old framework of testing for mediation, mediation was established by testing $b_{12}$ in equation 1, $b_{22}$ in equation 2, and $b_{32}$ in equation 3. So far, I have attempted to simulate values of $med$ and $Y$ that are consistent with values of the various regression coefficients using `rnorm` in `R`, such as the code below:               x <- rep(c(-.5, .5), 50)     med <- 4 + .7 * x + rnorm(100, sd = 1)           # Check the relationship between x and med     mod <- lm(med ~ x)     summary(mod)          y <- 2.5 + 0 * x + .4 * med + rnorm(100, sd = 1)          # Check the relationships between x, med, and y     mod <- lm(y ~ x + med)     summary(mod)          # Check the relationship between x and y -- not present     mod <- lm(y ~ x)     summary(mod)      However, it seems that sequentially generating $med$ and $Y$ using equations 2 and 3 is not enough, since I am left with no relationship between $X$ and $Y$ in regression equation 1 (which models a simple bivariate relationship between $X$ and $Y$) using this approach. This is important because one definition of the indirect (i.e., mediation) effect is $b_{12}-b_{32}$, as I describe above. Can anyone help me find a procedure in R to generate variables $X$, $med$, and $Y$ that satisfy constraints that I set using equations 1, 2, and 3?


title: Wigan scores after 30 minutes. Calculate home, away and draw in percentage terms with poisson regression

text: I have the following exercise: > Wigan v city 90 mins match: Wigan scores after 30 minutes. Calculate home, > away and draw in percentage terms with poisson regression, assuming pre- > match odds were 5.00 - 4.00 - 1.80 for Win/Draw/Win respectively. Please advise how to proceed with this question for my exam. ADDITION: after reading answers below and trying to answer this I got as far as below. Pls help. > (NOTE: I just kept it simple and assumed 2/3 of the game left) > > wigan mean, avg goals per game = 1.20 manu mean = 2.20 > > 60 mins left, so new wigan mean 2/3 x 1.20 = 0.8 new manu mean 2/3 x 2.2 = > 1.47 > > WIGAN WIN ODDS > > P(wigan scores 1 or more) = 1 - Possion(0, 0.8) = 1 - 0.45 = 0.55 > > OR, P(MANU FAIL TO SCORE) = Poisson(0, 1.47) = 0.23 > > WIGAN WIN ODDS = 0.55 + 0.23 = 0.78 > > Quote: WHAT ABOUT P(wigan 2 AND manu 1) = P(wigan goals > manu)??? > > MANU WIN ODDS > > Poss(manu 2 or more) = 1 - [Poisson(0) + Poisson(1)] > > = 1 - [0.23 + 0.34] > > = 0.43 > > DRAW ODDS > > X = ManU scores exactly 1 and Wigan scores exactly 0 from > > here on, AND manu scores same as wigan from here on. > > PossM(1) x PossW(0) = 0.34 x 0.45 = 0.153 > > P(manu = wigan in scoring from here on), 1.47 - 0.8 = 0.67 > > Poss(0 goals given mean = 0.67) = 0.51 > > So P(X) = [0.34 x 0.45] x 0.5 = 0.077 > > 0.78 + 0.43 + 0.077 = 1.28 > > 1.28 = 100% So, 0.78/1.28 x 100% = 61% > > So, 0.43/1.28 x 100% = 34% > > So, 0.077/1.28 x 100% = 5% > > DECIMAL ODDS = 100%/61% = 1.64 > > = 100%/34% = 2.94 > > = 100%/5% = 20.00 > > Clearly I made a mistake somewhere because the odds are too long on draw at > 19/1 even though draw should be a favorite or close to wigan win here.


title: Quantifying the relationship between two disparate time series

text: I have two time series that have a roughly similar trend, though both variables are noisy. This graph shows means and standard errors throughout a season of measurements. ![enter image description here](http://i.stack.imgur.com/7ThlT.png) I'd like to be able to make a quantitative statement about the relationship between these two data sets. While the two data sets were collected from the same experimental plots, the individual samples from which the means and standard errors were calculated are not meaningfully paired with one another, and you can see that the carbohydrate data set was measured more frequently. By taking a subset of the carbohydrate measurements that are closest to the microbial biomass measurement dates, I can make a scatterplot showing the means and standard errors that I think gives a fair visual representation of the relationship (TRS.ml is the carbohydrates): ![enter image description here](http://i.stack.imgur.com/kLvc0.png) This is where I am stuck. I'm not sure how to estimate regression coefficients or calculate an r2 value for a regression of this sort where I have estimates of uncertainty for both variables. Here are some approaches I have been considering:   1. Deming regression. I'm not sure that this would be the right approach. It seems to be more for data sets in which the same technique was used for both variables. If it is, my question is how would I calculate the variance ratio based on the information I have?   2. Regression of all underlying data points. This doesn't really work because the data are not meaningfully paired, so of the 80 or so microbial biomass measurements that underlie the data shown in the graphs here, I can't directly match them to individual measurements of carbohydrates. Matching them arbitrarily seems bad.   3. Regression of carbohydrate means by date against microbial biomass means by date. Basically regress the points in my scatterplot above but throw out the information about the uncertainty. This gives a high r2 driven by the coinciding peaks on July 1st, but to me, seems to overestimate the strength of the relationship.   4. Regression of all microbial biomass values against carbohydrate means by date or vice versa. This allows more of the underlying uncertainty to be incorporated while not forcing the pairing of unrelated data points in an arbitrary way. Again though, it does not incorporate the uncertainty in both variables. My question is which of these approaches, or any other unlisted approaches, would you recommend for quantifying the relationship between these two time series?


title: Free PDF for Bayes

text: Is there an good book/pdf similar to `Elements of Statistical Learning` that's available for free, online that deals with Bayesian statistics, ideally with code for `R`?


title: Why can't I approximate sinusoid with tanh activation function?

text: I've been playing with this applet for learning about backpropagation:   http://neuron.eng.wayne.edu/bpFunctionApprox/bpFunctionApprox.html What continues to confuse me, is why I can't approximate $sin()$ function with "Bipolar Sigmoid" activation function (supposedly $tanh$) while it works well with "Unipolar Sigmoid" (supposedly logistic function). Or is it just bug in the applet? Steps to reproduce:   1. Click on **Function** button and choose **Sin**.    2. Click on **Neural Network** and choose **Unipolar Sigmoid** for **Hidden Layer Activation Function**.   3. Click on **Training** and choose **1000** for **Number of Training Cycles**.   4. Click on **Train** few times to see it converge.   5. Then click on **Neural Network** and change activation function back to **Bipolar Sigmoid**.   6. Click on **Randomize Weights** to initialize training.   7. Click few times on **Train** to observe, that nothing much happens.


title: Performing binary logistic regression with equal number of cases and non-cases

text: The best way to ask my question is to present an example scenario: Let's say that the outcome of interest is lung cancer (1 = lung cancer; 0 = no lung cancer) and a researcher has 200k records (where 20k patients have lung cancer (cases) and 180k patients do NOT have lung cancer (non-cases)). Since only 10% of patients (20/200k) in the sample data have lung cancer, a researcher uses a random sample of 20k from the patients that do NOT have lung cancer. By doing so, the researcher would have a sample of 20k patients with lung cancer and 20k patients without lung cancer in their sample (the sample is reduced from 200k to 40k records). Are there any benefits to performing binary logistic regression with equal number of cases and non-cases when the actual distribution of the outcome is not equal? Or does this bias model estimates/predictive power? Thanks in advance!


title: Joint distribution of mid p-value and p-value

text: I have a question about the joint distribution of the mid p-value and p-value. We know that, for right tailed test with discrete test statistic $X$ with distribution $F$, the p-value is defined as $P=Pr(X \geq observed~X)$ and the mid p-value is defined as $mid~P=Pr(X \gt observed~X)+\frac{1}{2}*P(X=Observed)$. I would like to find out $Pr(P_{mid} \leq t~and~P \gt t)$ for some $t \in (0,1)$. Here $P$ is the p-value and $P_{mid}$ is the mid p-value.


title: Elementary statistics for jurors

text: I have been summoned for jury duty. I am conscious of the relevance of statistics to some jury trials. For example, the concept of "base rate" and its application to probability calculations is sometimes - perhaps always - relevant. What statistical topics might a person in my situation usefully study, and what materials would be suitable for somebody with my background? I have a "hard science" degree and so have some limited statistical knowledge, but my skills are rusty. I work full time and don't have a lot of time before my jury duty. So it would be appropriate for the answers to focus on elementary concepts, simple problem solving skills and their application to relevant problems (and the limitations of those concepts and methods, of course).


title: Forecasting with arimax model including xtransf

text: I am trying to forecast with an arimax model including the `xtransf` argument. I used the example given in chapter 11 of the book Time Series Analysis With Applications in R, Jonathan D. Cryer & Kung-Sik Chan, 2nd ed. So this is basically an intervention Model. I just did not include the `xreg` argument, so the Additive outliers are not included for the moment. Here are the codes and the Error I got.               >library(TSA)     >data(airmiles)      >air.m1=arimax(log(airmiles),order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12),xtransf=data.frame(I911=1*(seq(airmiles)==69),     I911=1*(seq(airmiles)==69)),transfer=list(c(0,0),c(1,0)),method="ML")     > air.m1          Call:     arimax(x = log(airmiles), order = c(0, 1, 1), seasonal = list(order = c(0, 1,          1), period = 12), method = "ML", xtransf = data.frame(I911 = 1 * (seq(airmiles) ==          69), I911 = 1 * (seq(airmiles) == 69)), transfer = list(c(0, 0), c(1, 0)))          Coefficients:               ma1     sma1  I911-MA0  I911.1-AR1  I911.1-MA0           -0.5379  -0.7644   -0.1290      0.8901     -0.2419     s.e.   0.0854   0.1532    0.0606      0.1239      0.0513          sigma^2 estimated as 0.0009739:  log likelihood = 199.42,  aic = -388.84          >predict(air.m1,n.ahead=20)     Error in dim(data) <- dim : attempt to set an attribute on NULL      Any idea where my mistake is?


title: Hypothesis testing using the likelihood ratio test

text: How can I apply the likelihood ratio test when I have the following hypothesis for collocation discovery: $$H_{1}: P(w^{2}|w^{1}) = p = P(w^{2}| \neg w^{1})$$ $$H_{2}: P(w^{2}|w^{1}) = p_1 \neq p_{2} = P(w^{2}| \neg w^{1})$$ Basically, $H_{1}$ describes the hypothesis that the second word ($w^{2}$) in a possible collocation is independent from the first one ($w^{1}$) and $H_{2}$ is the opposite, that the second word depends on the first one. It's assumed that these bigrams follow a binomial distribution $$b(k;n,x) = \displaystyle \binom{n}{k}x^{k}(1-x)^{(n-k)}$$ However, the book I'm following says that the likelihood for $H_{1}$ is $$L(H_{1})=b(c_{12};c_{1},p)b(c_{2}-c_{12};N-c_{1},p)$$ and for $H_{2}$ is $$L(H_{2})=b(c_{12};c_{1},p_{1})b(c_{2}-c_{12};N-c_{1},p_{2})$$ where $c_{1}$, $c_{2}$ and $c_{12}$ is the number of occurrences of $w^{1}$, $w^{2}$ and $w^{1}w^{2}$, respectively. Why is that the case? I understand that the likelihood is defined as $L(p, x)$ where $p$ is a given parameter and $x$ is some observation or data, which is the opposite question one would ask when dealing with a probability mass function $f(x, p)$, that is, the probability of having $x$ given some parameter $p$. Moreover, I also know that when we have a probability mass function of a series of observations $f(x_{i}, p)$ then, assuming that the observations are independent from each other, the likelihood of a parameter $p$ given those observations is $L(p, x_{i}) = \prod_{i} f(x_{i}, p) $, but I don't comprehend why the example described above apparently is multiplying the probability of getting a $w^{1}w^{2}$ and $\neg w^{1}w^{2}$ in both hypothesis. Can someone explain that to me?


title: Is is logical to have AR terms in stationary series?

text: Ok so the theory says that ARMA should be applied to stationary series only. By definition stationary process is one with constant mean and variance. However if we use AR terms to model the series wouldn't that lead to dependence of the the series on past value? Which would in turn lead to non- constant mean and variance? Also related question that has not yet been answered: ADF Test on Volatile Series. Expected results?


title: linear regression with 3-dimensional points

text: I've been hunting around looking for a way to do this problem. I need to create a program to calculate linear regression for 100 3-dimensional points. I also have the matching outcomes of the points, so it's like a training set rather than a testing set. I'm also instructed to consider the bias term, but I'm not sure what that means. I can find good documentation for 2-dimensional points, but not for 3 dimensions. Since there's so many equations for doing XY data points, surely there's a similar equation set for XYZ points. I'm not really desiring to be given the equation outright, but more interested in an understanding of the function and how it works, as well as how it's determined. Thanks to anyone who read and who can help me with this problem.


title: Breiman's RF code

text: Is anyone familiar with usage of Random forest code Breiman and Cutler (http://stat-www.berkeley.edu/users/breiman/RandomForests/cc_home.htm).?? I wanted to know how to use that FORTRAN code to predict the class labels of records in a unlabeled dataset?


title: ARIMA with seasonality in Statsmodels

text: I'd like to have a _seasonal_ ARIMA model implemented with Statsmodels ARIMA. Specifically, I'd like to log before the weekly seasonality and then be able to make forecasts. Perhaps an example with ARIMA's from_formula method could accomplish this. I'd also love to be able to do this with patsy. Here's my sample code for logging before the weekly seasonality, and then transforming back to compare to the original time series (I've also skipped checking the validity of the model through testing stationarity and the residuals):               import pandas as pd     import numpy as np     from statsmodels.tsa.arima_model import ARIMA          # ts is a time series     logged_ts = np.log(ts)     # Differencing by the week forces us to drop the first 7 values.     diffed_logged_ts = (logged_ts - logged_ts.shift(7))[7:]          p = 0     d = 1     q = 1          arima = ARIMA(diffed_logged_ts.values, [p, d, q], exog=None, dates=diffed_logged_ts.index, freq='D', missing='none')     diffed_logged_results = arima.fit(trend='c', disp=False)     predicted_diffed_logged = diffed_logged_results.predict(exog=None, dynamic=False)     predicted_diffed_logged_ts = pd.Series(predicted_diffed_logged, index=diffed_logged_ts.index[d:])     predicted_diffed_logged_ts = np.exp(logged_ts.shift(7) + diffed_logged_ts.shift(d) + predicted_diffed_logged_ts)          concatenated = pd.concat([ts, predicted_diffed_logged_ts], axis=1, keys=['original', 'predicted'])     print concatenated[-7:]      What do you think of this approach? I hope there's a less error-prone way coming in a future version of Statsmodels. Could someone tag this question with "statsmodels"? Thanks!


title: Random slopes and intercepts using glmmadmb in R - negative binomial

text: I am new to mixed modeling and using the glmmadmb package. I am interested in using glmmadmb to model count data with a negative binomial distribution, using year and site sampled as my random effect variables, as both should be allowed to vary independently. I would like to consider random intercepts for both factors, but also random slopes for sites across all years. I am not sure how to code this in R, and have tried several different ways. Usually R just takes a super long time and does not seem to complete the computation - which leads me to think my syntax is incorrect. Any suggestions/ thoughts would be greatly appreciated! I'm using Windows 7 with RStudio v.0.98.312. Here is what works - random intercepts model:                  mod <- glmmadmb(species~(1|year)+(1|site),data=cs,           family="nbinom2",link="log")      Here are some iterations I have tried which did not seem to work:                 mod <- glmmadmb(species~(1+site|year)+(1|site),data=cs,...)       mod <- glmmadmb(species~(site|year)+(1|site),data=cs,...)       mod <- glmmadmb(species~(1|year)+(1|site)+(0|site),data=cs,...)      Thank you, Tiffany


title: Data Correlation for Predictive Basketball Model

text: **Background Information: (scroll for question)** Attempting to find correlation between the average number of fouls committed by each basketball team in a game and the referees calling the game. Such a relationship would provide a predictive model for the total expected number of fouls predicted in any match. With the extraordinarily human dependent element of this data and the relative unpredictability of a basketball game, such a correlation would be weak at best, but interesting nonetheless. * * * My initial strategy was to consider what portion of the actual outcome came from the expectations of the teams and refs respectively, leading to the following relationship: $$aX + bY = c \\\ a = sum \ average \ commited \ fouls \ by \ both \ teams \\\ b = average \ num \ fouls \ called \ by \ specific \ refs \\\ c = actual \ outcome = total \ num \ fouls$$ I assumed $X$ & $Y$ would be simple probabilities summing to 1 to complete a total outcome. However I quickly realized this wasn't a possibility as any outcome $c$ greater than $a$ or $b$ would be invalid and thus disregarded in the predictive model, destroying accuracy in the process. * * * I have compiled data from 20 games in the google doc below. As expected for such a small sample size the data is anything but conclusive. That being said, data is available for 1000s of samples once a statistical model is selected. I've included a graph of the data below for convenience. ![sample data](http://i.stack.imgur.com/P8IWD.jpg) Sample Data: http://goo.gl/SrlRnQ [google doc] * * * **Question & Inquiry:** Why I've come to Cross Validated is to seek out an ideal statistical method for determine such a correlation in this data. It seems a liner regression may be the route to go although I'm not sure how to proceed. The available data is on the order of 1000s of sets, not just the 20 I have provided. Given this data, what is the best way to determine a consistent, or at least average, prediction of the contributions of both the teams and the refs in the outcome. I.e. what portion, on average, of the outcome can be directly attributed to each party across 1000s of samples. My initial model led me to believe an average, universal scaler could be calculated for both a and b given above to predict the outcome. As this would be a mere average, it wouldn't be terribly accurate, but I'm curious to see just that the correlation is.


title: How is adjusted coefficient of determination ($R^2$) linked to the F-values of a test against zero when adding a new variable?

text: Someone claims that the adjusted $R^2$ will increase with the addition of an extra variable. I wonder why, as it is called adjusted (in contrast to the normal $R^2$). The only condition it has to satisfy (to increase the adjusted $R^2$) is that the F-value (by the way, how is it simple to calculate it?) of the null hypothesis that the new variable is greater than 1. Can someone give me a hint where the links between the adjusted $R^2$ and the F-Stat of that test are? And however, who wants to include a new variable in a multiple OLS regression model anyway, if the beta was tested to be 0? Therefore adjusted $R^2$ always changes.


title: Interpret the pairwise comparison plot from SPSS

text: For my user study, I am trying to find out to what extent people prefer to move by foot. They had the option to rank their preference for 'walking', 'jogging' and 'running' for a several amount of questions. Now I want to interpret the data by means of a Kruskal-Wallis test (`Analyze -> Non Parametric Test -> Independent samples`). The following analysis gives me pairwise comparisons. However, I do not understand how the pairwise comparisons' plot has been created when I look at the table below it. How is the plot from the pairwise comparison generated when looking at the table below the plot? ![table overview](http://i.stack.imgur.com/mtIqL.png) ![rejection of null hypothesis](http://i.stack.imgur.com/WeNUy.png) ![boxplots](http://i.stack.imgur.com/l1dUj.png) ![Pairwised comparison plot](http://i.stack.imgur.com/EobhW.png)


title: exploratory analysis of associations between many variables

text: Can anyone suggest good references or examples of testing for associations between many (50-100) different variables/features (which, in my case, are primarily categorical or binary) This seems to me to be related to a) association rule / frequent itemset mining b) bayesian network learning The general idea is to systematically make lots of 2x2 contingency tables and run pairwise tests of association (adjusting for multiple testing). This could be extended to all the 2-way tests conditional on each of the other variables, and to 3-way tables, etc. However, I'm wondering if there's been (but have been surprisingly unable to find much) prior work on: * what are the best metrics to use for assessing associations in this context  * software for efficiently exploring and visualizing these spaces of associations  * efficient implementations for running all these tests  * properly controlling for multiple testing, when some features are correlated (whether known or not)


title: Random variables: True or False

text: I sometimes get confused when reading statistical definitions when they mention random variables (RV). Are they talking about a single draw? Are they talking about an estimator? To make things clear for myself (and others reading the question) I have created a list of statements that I think are true. Is there anything that I have said here that is false? 1) An element in a population **IS NOT** a RV since it cannot take on a set of different values. 2) A parameter of a population (e.g. mean) **IS NOT** a RV since it is based on elements which are not RVs. 3) A random draw from a population **IS** a RV since its realization is based on the proportion of that element in the population. 4) A sample is made of RVs since each element of the sample is itself a RV. 5) A statistic of a sample (e.g. sample mean) **IS** a RV since the sample is made up of RVs. 6) An estimator of a population parameter (e.g. sample mean) **IS** a RV since it is based on a sample and the sample is made up of RVs. 7) A sequence of statistics (e.g. sequence of sample means) is a sequence of RVs. 8) A sequence may converge to either to a constant or a RV.


title: How do you correlate binary & ordinal variables?

text: Which test to use to calculate correlations between them?   * Cramer's V,   * Kruskal-Wallis,   * or something else?


title: How to test if time series are not cointegrated?

text: Is there a test whose null hypothesis is that time series are cointegrated? There are lots of tests whose null is no cointegration. How can I calculate a p-value for the null that two series are cointegrated? I assume a p-value of 0.95 for a null of no-cointegration is not the same as a p-value of 0.05 for a null of cointegration.


title: Stationarity consideration in ARIMA using KPSS test

text: I have data, which I am sure has a downward trend. I am trying to forecast this data using ARIMA and I want ARIMA to consider the trend when it is forecasting. The first step in ARIMA is to determine the order of differencing for which I am using the KPSS test. The KPSS test in R has an option 'trend' and 'level'. Which one should I choose if I know that my data has a downward trend? I am confused because I dont understand what is the input ARIMA is expecting if trend has to be taken into consideration for forecasting. Another question: Is there any sure shot way of knowing if the data follows a trend?


title: Raw return vs. percentage return to calculate volatility

text: I am using squared return as a proxy to calculate volatility, however i'm not sure whether to use raw return or percentage return. Under raw return all return estimates are below 1, however under percentage return there is a mix of return greater than 1 and less than 1. Percentage return below 1 would end up as a volatility figure less than the percentage return itself, on the other hand percentage return above 1 would end up as a volatility figure greater than percentage return. My question is: Doesn't this pose a problem when calculating volatility in that there is an over estimation when the return is above 1? I am going to use the data to fit an ARMA-GARCH model, would there be any difference if I used percentage or absolute values?


title: Satterthwaite vs Kenward-Roger approximations for the df in mixed effects models

text: The `lmerTest` package provides an ANOVA function for linear mixed effects models with optionally Satterthwaite's (default) or Kenward-Roger's approximation of the degrees of freedom. What is the difference between these two approaches? When to chose which?


title: User intent identification in web search queries using unlabeled data

text: Given a list of unlabeled search queries, is it possible to say which query is related to a certain topic? For example, if we want to select the queries where the intent of the user is to download or listen to music online, and the dataset contains queries like these:   * best songs in 2013   * facebook   * listen to music online   * download free music   * youtube   * news   * college applications   * calculus tutorial   * vampire diaries I found a few papers but they use supervised or semi-supervised methods or extra information such as click counts. Thanks in advance.


title: Regression with Census variables

text: I seek to model change in school performance over three years (2010 to 2012) using a number of socio-economic variables. The socio-economic data is available from two census periods: 2000 and 2010. I would like advice on what option is better for me; regressing the performance for each of these years against the socio-economic data of 2010, or regressing 2011 and 2012 performance with estimated census data using the change between 2000 and 2010 census data (i.e. only 2010 performance uses 2010 data).


title: Missing data in repeated measures design

text: I have data from a a simple _2(voice system, manual system)_ x _2(easy, hard)_ within subjects experiment. Multiple DVs were collected. Due to issues with both the systems tested, data for one system condition or another is corrupt (so badly that it can be considered missing) for many subjects (ss). Out of data for 34ss, only 14 have data for both system conditions. The remaining 20 have data for the voice system (9) or the manual system (11). To list it out:   * 14 x data for both    * 11 x data for manual system condition, but missing voice   * 09 x data for voice system condition, but missing manual **I am considering my options:**   * I can use a repeated measures analysis with listwise deletion of missing data and run my analysis on only the 14 with full data.   * I can use multiple imputation to attempt to fill in for 20 missing sets of data, however, that is a lot of data to impute. I'm wondering if that is wise.   * I can redesign the experiment into a between/within ss experiment, randomly discarding data from the 14 who have both systems dara and arriving at a balanced 22n in each newly coined 'between' group. The difficulty group would still be analyzed as within.   * I can begin collecting again. Sadly, in repairing the facility it was also 'upgraded', meaning enough changes that I'll have to start from scratch. Any opinions on which is the best path forward? Any paths I have not considered?


title: What is linear separability of classes and how to determine

text: This question may seem too trivial but my basics are not strong and I shall appreciate help in these concepts. For an n dimensional feature vector and 3 class problem does linear separability need to be checked?,   1. I am familiar with XOR problem which cannot be modeled by neural network since the class is not linearly separable. What exactly do we mean by linearly separable? If the class boundary can be drawn as a curve and we can divide the patterns into 2 classes then can it be called linearly separable?   2. How to decide if the classes are linearly separable?    3. Does linear separability apply to more than 2 class problem? Say for 3 classes, how does one decide? Thank you.


title: How to create item parcels in SPSS Amos?

text: I am trying to run a SEM to see the causal relationship among latent variables. I ran the CFA and SEM using Amos, but I realized that the model is very complicated since one factor, which I originally thought to be unidimensional, is two-dimensional, and that makes the model look very messy. Since that factor is not derived from a theoretical background, it seems hard to force to use it as unidimensional. Also, a construct validity test also indicated that the factor is two-dimensional. I know that creating parcels and running the path model has lots of limitations, such as not including error, and it is not necessarily used to create a parsimonious model, but I want to compare whether there are any big differences between CFA and a model with item parceling. I saw lots of articles that have debated item parceling, but couldn't find how to conduct item parceling in Amos. How can I conduct item parceling in Amos? Thanks!


title: Interpreting coefficients for random effects models with extremely unbalanced data

text: I'm currently working with a data set that has numerous samples collected over time at different sites in a study area, and I'm interested in detecting a trend over time for that area. I know that in an ideal experimental or balanced situation, using a random slope and intercept model is a great way to get at the overall trend within the study area. With our data, however, many of the sites are missing samples and a handful of the sites only have one data point. I'm curious if there's a way to intuitively understand how the sample imbalance will affect the estimate of the overall slope? To put it differently, are there ways to know if sample imbalances are causing problems, or are there things I can look for in my model output that would indicate I shouldn't trust what the model is estimating? I created a contrived example with 20 data points to look at this. I put 10 data points with a slope of 1 into one site (a), and put the other 10 data points with a slope of -1 into unique sites (b through l). I had assumed that when I looked at both a random intercept and random slope and intercept model that they would be somewhat similar, or that at least the latter would give more weight to the site with good data over time.               > library(lme4)     > set.seed(9999)          > x = c(0,1,2,3,4,5,6,7,8,9,0,1,2,3,4,5,6,7,8,9) + rnorm(20,mean=0,sd=0.1)     > y = c(0,1,2,3,4,5,6,7,8,9,9,8,7,6,5,4,3,2,1,0) + rnorm(20,mean=0,sd=0.1)     > z = c(rep('a',10),'b','c','d','e','f','h','i','j','k','l')     > z = factor(z)          > m0 = lm(y~x)     > m1 = lmer(y~x+(1|z))     > m2 = lmer(y~x+(1+x|z))          > summary(m0)     > summary(m1)     > summary(m2)     > anova(m1,m2)      As expected, the slope of the linear model was near zero, but the results for the two mixed effects models were nearly opposite. Even though sites b through l only have one data point, it seems like they contribute more towards the slope because the trend is occurring over so many sites. The random slope and intercept model was also preferred to using model selection criteria.                > summary(m0)$coefficients                     Estimate Std. Error    t value    Pr(>|t|)      (Intercept)  4.53784796  1.2586990  3.6051890 0.002023703      x           -0.01178748  0.2335094 -0.0504797 0.960296079           > summary(m1)      Linear mixed model fit by REML ['lmerMod']      Formula: y ~ x + (1 | z)            REML criterion at convergence: 62.0877            Random effects:       Groups   Name        Variance Std.Dev.       z        (Intercept) 33.30788 5.7713         Residual              0.01583 0.1258        Number of obs: 20, groups: z, 11           Fixed effects:                  Estimate Std. Error t value      (Intercept) -0.03597    1.74163   -0.02      x            0.99332    0.01386   71.66           Correlation of Fixed Effects:        (Intr)      x -0.036           > summary(m2)      Linear mixed model fit by REML ['lmerMod']      Formula: y ~ x + (1 + x | z)            REML criterion at convergence: 31.0386            Random effects:       Groups   Name        Variance Std.Dev. Corr        z        (Intercept) 7.78818  2.7907                   x           0.37691  0.6139   -1.00       Residual             0.01524  0.1234              Number of obs: 20, groups: z, 11           Fixed effects:                  Estimate Std. Error t value      (Intercept)   8.2121     0.8566   9.587      x            -0.8201     0.1882  -4.358           Correlation of Fixed Effects:        (Intr)      x -0.999           > anova(m1,m2)      Data:       Models:      m1: y ~ x + (1 | z)      m2: y ~ x + (1 + x | z)         Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)          m1  4 66.206 70.189 -29.103   58.206                                   m2  6 36.745 42.719 -12.372   24.745 33.462      2  5.419e-08 ***      ---      Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      I see that under this extreme example, the random slope and intercept have an almost perfect correlation. Is what I can pull from this is that, in a sense, the model gives more value to the sites with only one data point because the overall trend is so strong but across multiple sites, but that I should view the slope estimate this model produces as suspect with such a high correlation? Is there anything else that should look for? For my specific study, I could also set some sort of criteria for what level of replication I thought was necessary to make proper inferences, e.g. eliminate all the sites that less than five samples. Many thanks for your thoughts.


title: Difference between GLS and SUR

text: I've been reading some about Generalized Least Squares (GLS) and trying to tie it back to my basic econometric background. I recall in grad school using Seemingly Unrelated Regression (SUR) which seems somewhat similar to GLS. One paper I stumbled on even referred to SUR as "special case" of GLS. But I still can't wrap my brain around the similarities and differences. so the question: What are the similarities and differences between GLS and SUR? What are the hallmarks of a problem which should use one method over the other?


title: Feature relationship based class separability

text: I am a computer science guy, not a mathematician so kindly excuse me if there is any ridiculous error in my problem description. I have two clusters $C_1$ and $C_2$ in a feature space spanned by $k$ features scaled between $0$ and $1$. That means each feature ranges between $0$ and $1$ after min-max normalization of the feature values. The clusters are generated while sampling from two separate multivariate Gaussian distributions $N_1(f_1, f_2, ...f_k )$ and $N_2(f_1, f_2, ...f_n )$ where $f_i \in [0,1]$, over the scaled features. I do not have a condition on means or covariance matrices. However it is given that any pair of vectors (or instances) within a given cluster should not have the Euclidean distance more than a specified scaler $D$. Now two arbitrary vectors $A$ and $B$ are picked up from clusters $C_1$ and $C_2$ respectively. Now $R_A$ and $R_B$ are the feature orderings based on their corresponding values in the vectors (for example, if $A$=(0.8,0.2,0.7), $R_A$ is an ordered list $(f_2<f_3<f_1)$). What is the probability that _Kendall's tau_ distance between $R_A$ and $R_B$ is $0$. I need an expression for $Pr(K(R_A,R_B)=0 | D=d)$. I hope I could explain the problem unambiguously. However, kindly let me know if there is any confusion. Thanks.


title: Understanding the mixed effect model thru creating my own optimization-based model

text: Following Macro's suggestion, I am creating a new thread for my questions. (It's somehow related to the original question, which had been discussed here: Understanding the linear mixed effects model equation and fitting a random effects model with weights in R Specifically, I am trying to understand a random intercept model in the form of $${\bf y}={\bf X}{\boldsymbol \beta} + {\bf Z} {\bf b} + {\boldsymbol \varepsilon} $$ (Following Macro's notation for ${\boldsymbol y_i}$ in the previous question, I stack the group formulas together vertically and obtained one big formula for the whole dataset with all groups). My questions are:   1. In a random intercept model (all group sharing the same beta), will the following situation arise: all data points in a group reside on one side (the same side) of the regression line of that group? Here is my thinking: even though for a certain group, the overall beta is a very crappy one for that group. By adjusting the intercept for that group, I at least still can get the regression line passing the cloud, am I right? Here by "passing the cloud" I mean not all points will reside on the same side of the regression line.   2. In a random intercept model, for each group, which point does the regression line pass thru? In the general OLS, we knew that it's the point $(\overline{x}, \overline{y})$ that is passed thru by the regression line. But for each group in a random intercept model, which point is that "central" point?    3. Moreover, in `R`, what's a convenient way to visualize what's happening within a certain group in a mixed model?    4. How do I understand the ${\bf b}$'s? What are these ${\bf b}$'s? May I say they are simply the $\overline{y}_i - \overline{y}$, i.e. the vertical axis of the group centers - the vertical axis of the whole-data center?   5. Why does `lme` explicitly impose that the sum of ${\bf b}$ is $0$ ?   6. To understand the `lme` models better, I am creating my own toy models and doing some comparative studies: To cook my own model, I decide to solve the following least-square problem: $$ \min \left( \sum_{k=1}^{K} w_{1k} {\bf b}_k^2 + \sum_{i=1}^{N} w_{2i} \varepsilon_i^2 \right) $$ here $\\{ {\bf b}_k \\}_{k=1}^{K}$, where $K$ is the number of groups. $\\{ \varepsilon_i \\}_{i=1}^{N}$, where $N$ is the number of individuals. $w_1$ is a set of weights on the individuals and $w_2$ is a set of weights on the groups at the group level. My purpose is to create a model that resembles the random intercept model, but ditching the Gaussian distribution assumption for the ${\bf b}$'s and I would like to have the regression line in each group(sharing the same beta) passing the data cloud for that group (Please ref. Q1 above).   * Do you think my model makes sense?    * Could you please help me critique my model vs. the `lme` random intercept model?    * Shall I add more constraints such as $\sum {\bf b}_k = 0$and $\sum \varepsilon_i =0$? Thanks a lot for your help!


title: How to produce a pretty plot of the results of k-means cluster analysis?

text: I'm using R to do K-means clustering. I'm using 14 variables to run K-means   * What is a pretty way to plot the results of K-means?    * Are there any existing implementations?   * Does having 14 variables complicate plotting the results?  I found something called GGcluster which looks cool but it is still in development. I also read something about sammon mapping, but didn't understand it very well. Would this be a good option?


title: How to estimate how many people attended an event (say, a political rally)?

text: A student asked me today, "How do they know how many people attended a large group event, for example, the Stewart/Colbert 'Rally to Restore Sanity' in Washington D.C.?" News outlets report estimates in the tens of thousands, but what methods are used to get those estimates, and how reliable are they? One article apparently based their estimate on parking permits... but what other techniques do we have? Please note I am not talking about capture/recapture experiments or anything of the like. I don't have any idea. I would guess in advance that there aren't specific methods for something like this, and whatever's there is very ad hoc (such as how many parking permits were sold). Is this true? For purposes of national security - of course - it would be possible to have an analyst sit down with satellite photographs and physically bean-count the number of people there. I doubt this method is used very often.


title: Question about KNN algorithm

text: 1. Does KNN algorithm have a centroid as _k_ -means?   2. Is there a way to obtain the centroid for the classified data by KNN?   3. Is there a way to compare SVM classification with KNN classification?


title: Are inconsistent estimators ever preferable?

text: Consistency is obviously a natural and important property estimators, but are there situations where it may be better to use an inconsistent estimator rather than a consistent one? More specifically, are there examples of an inconsistent estimator which outperforms a reasonable consistent estimator for all finite $n$ (with respect to some suitable loss function)?


title: Anderson-Darling vs Kolmogorov-Smirnov (Many points in sample)

text: ### Previous question In my previous question (Kolmogorov-Smirnov test - reliability) I misused two- sample KS test for normality testing of one-sample data. I was given an advice that I should use one-sample KS test, also it was recommended to me to use Anderson-Darling test. ### Current situation I use following C++ implementations of Anderson-Darling (AD) and Kolmogorov- Smirnov (KS) tests to check if given cluster has normal distribution or not. ### Problems If i deal with lets say hundreds of points per cluster, AD gives me satisfactory results. Nevertheless, If I have thousands points in my cluster, I receive large values of A-squared from AD that are far away from AD's critical values. I have read here that AD fails for large number of points and that it is better to use KS. However, I am not sure, if I did something wrong, but KS is equal to 1.0 in almost every case (cross checked in Matlab too). ### Example data Example data are located here. It represents a cluster of gaussian mixture points. I want to test this cluster for normality (and I want to fail in that test for this data). Consequently, I want to split data to sub-clusters, test each sub-cluster again and split again until I end up with "leaf clusters" that are normally distributed.


title: Why can't one generalize the Kolmogorov-Smirnov test to 2 or more dimensions?

text: The question says it all. I've read both that one can't generalize KS to a dimension equal or larger than two, and that famous implementations like that in Numerical Recipes are simply wrong. Could you please explain why is so?


title: Providing the actual measured error instead of weights to the lm function in R (or similar linear fit function)

text: Is there any way to use the linear fit function or R with my own errors instead of the value weights? That is, I have a physical experiment in which we know the errors for each point, for instance:               > df         x        y   yerror     1   0 3097.460 55.65483     2  60 4624.340 68.00250     3 121 7013.675 83.74769          > dput(df)     structure(list(x = c(0, 60, 121), y = c(3097.46, 4624.34, 7013.675     ), yerror = c(55.6548290806827, 68.0024999540458, 83.7476865352112     )), .Names = c("x", "y", "yerror"), row.names = c(NA, -3L), class = "data.frame")      Is there any way to provide this to the **lm** function in R (or a similar one) but to give the actual errors instead of having R calculating the statistical ones? For instance, something like the **LINFIT** function from **IDL language** : http://www.astro.virginia.edu/class/oconnell/astr511/idl_5.1_html/idlea.htm I see I can give the weights as a parameter to **lm** , but what I want is to give the errors we measured. EDIT: I'm putting the code and the output in both R and IDL to see the differences. **IDL/GDL**               X = [0, 60, 121]     Y = [3097.46, 4624.34, 7013.675]     YERROR = [55.6548290806827, 68.0024999540458, 83.7476865352112]          result = LINFIT(X, Y, CHISQR=chi2, COVAR=varco, MEASURE_ERRORS=YERROR, SIGMA=sig)     print, result           3006.64      31.4443     print, chi2           24.1395     print, varco           2755.76     -29.0644          -29.0644      0.65546     print, sig           52.4954      0.80960      And now in **R** :               > lfit = lm(y~x, data=df, weights=yerror^-2)     > summary(lfit)          Call:     lm(formula = y ~ x, data = df, weights = yerror^-2)          Weighted Residuals:          1      2      3       1.632 -3.955  2.415           Coefficients:                 Estimate Std. Error t value Pr(>|t|)       (Intercept) 3006.640    257.920  11.657   0.0545 .     x             31.444      3.978   7.905   0.0801 .     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 4.913 on 1 degrees of freedom     Multiple R-squared: 0.9842, Adjusted R-squared: 0.9685      F-statistic: 62.49 on 1 and 1 DF,  p-value: 0.08011           > vcov(lfit)                 (Intercept)          x     (Intercept)  66522.7761 -701.59982     x             -701.5998   15.82237          > vcov(lfit)/4.913^2                 (Intercept)          x     (Intercept)  2755.98492 -29.066714     x             -29.06671   0.655508      As you can see the covariance matrix from the lm function in R, is multiplied by the square of the "Residual standard error". If that happens all the time, I guess I can just divide it by it and then calculate the error from the right covariance matrix. Not only that, but the presented error by R seems multiplied by the "Residual standard error" (almost, it's not really exact).


title: Creating a composite score out of a set of categorical variables

text: I am running an experiment that has a set of categorical variables (yes/no responses). If I wanted to reduce these categorical variables into a single score (e.g. a summed score/average score).   1. Is there an appropriate method for doing this?    2. I know that rasch/IRT models can be used to evaluate the scale properties. Can the rasch model be used as a basis for creating a composite score?    3. If so how?


title: Is machine learning useful for comparing test group and control group

text: If you do a hypothesis to test the effectiveness of a treatment, or a marketing campaign, you want to be sure that the two groups are comparable. You can compare some relevant quantities between the two groups. Even if these do not differ significantly, there may still be hidden differences between the two groups. Suppose we compare costumer spending. The test group are costumers whose names start with the letter a, b, c, d, or e. The time and amount at check out are stored. Can machine learning play a role in finding (statistically significant) differences between the groups? Are there any references on this kind of use of machine learning?


title: Is p-value involved in Neyman-Pearson Testing procedure?

text: It seems to me that, after having a predetermined level of significance, the Neyman-Pearson approach only requires to see if the observed x lies in the rejection region or not. In this case, how to relate it to p-value? Or is it an alternative way to judge if x is rejected or not by using p-values? After reading some references about NP and Fisher's treatment, I am confused. Thank you for any suggestion.


title: Growing Self-Organizing Map for mixed-type data

text: I am trying to write code to build a growing SOM for mixed-type data. I came across a paper Growing Self-Organizing Map with cross insert for mixed-type data (http://www.sciencedirect.com/science/article/pii/S1568494612001731). It was very interesting and handled both categorical and numeric data in a unified way. However, my dataset has variables/attributes that can have multiple values (for example: attribute "interests" can have more than one value - movies, sports, and so on...). I got stuck at handling such attributes. Any inputs how to handle attributes with sets of values in a mixed-type dataset? References to material that talk about this problem would be greatly appreciated.