title: General rules for forming groups for analyses

text: I am analysing EEG data. I have data from 19 electrodes from 25 participants divided into two groups (11 and 14 subjects). I asked my student to compare these two groups with some statistical test. Surprisingly, he calculated means for every electrode (giving two groups with 19 numbers in each). I suppose that it is more valid to calculate means for each participant (resulting in two groups with 11 and 14 numbers). Who is right? Or are we both right? What are the general rules for forming groups for analyses?


title: Moments of a distribution - any use for partial or higher moments?

text: It is usual to use second, third and fourth moments of a distribution to describe certain properties. Do partial moments or moments higher than the fourth describe any useful properties of a distribution?


title: Intuition / interpretation of a distribution of eigenvalues of a correlation matrix?

text: What is your intuition / interpretation of a distribution of eigenvalues of a correlation matrix? I tend to hear that usually 3 largest eigenvalues are the most important, while those close to zero are noise. Also, I've seen a few research papers investigating how naturally occuring eigenvalue distributions differ from those calculated from random correlation matrices (again, distinguising noise from signal). Please feel free to elaborate on your insights.


title: How to identify points and an unknown distribution in a two type clustering problem?

text: I have a data set that contains two types of points. The first type of points come from an N(0,1) distribution. The second type of points come from an N(m,v) distribution for some real m and some positive, real v. The objective is to classify each point as type 1 or type 2, and to identify m & v. We have no apriori information about m & v. Any ideas?


title: What is the distribution of OR (odds ratio)?

text: I have a bunch of articles presenting "OR" with a- 95% CI (confidence intervals). I want to estimate from the articles the P value for the observed OR. For that, I need an assumption regarding the OR distribution. What distribution can I safely assume/use?


title: OpenCV parameters of Random Trees

text: OpenCV provides an implementation of random forest named random trees and derived from a decision tree class. One parameter to train the random forest is the maximum depth, which in the provided examples is typically between 10 and 20. I learned that random forest is generally grown to its full depth and no pruning is done and, therefor, other random forest implementations do not provide a parameter of maximum depth. Why then, without consulting the source code, which might give the answer, does the OpenCV implementation provide this parameter and is it meaningful to limit the maximum depth in a random forest?


title: What types of situations would linear regression outperform randomforest?

text: I'm working on a data set and linear regression is beating a randomforest (the metric is weighted mean absolute error - I have a list of weights and each difference $|\text{predicted} - \text{actual}|$ is multiplied by that rows weight (which is a $\sqrt{\text{time}}$ so longer time between observations are charged a higher penalty)). I find this a bit surprising in light of past experience. I was curious if there are types of data that this would typically happen with AND if there are any types of algo/ml technique that's even better than linear regression.


title: Assessing forecastability of time series

text: Suppose i have a little over 20.000 monthly time series spanning from Jan'05 to Dec'11. Each of these representing global sales data for a different product. What if, instead of computing forecasts for each and every one of them, I wanted to focus only on a small number of products that "actually matter"? I could rank those products by total annual revenue and trim down the list using classical Pareto. Still it seems to me that, although they do not contribute much to bottom line, some products are so easy to forecast that leaving them out would be bad judjement. A product that sold 50$ worth each month for the past 10 years might not sound like much, but it requires so little effort to generate predictions about future sales that I might as well do it. So let's say I divide my products in four categories: high revenue/easy to forecast - low revenue/easy to forecast - high revenue/hard to forecast - low revenue/hard to forecast. I think it would be reasonable to leave behind only those time series belonging to the fourth group. But how exactly can I evaluate "forecastability"? Coefficient of variation seems like a good starting point (I also remember seeing some paper about it a while ago). But what if my time series exhibit seasonality/level shifts/calendar effects/strong trends? I would imagine I should base my evaluation only on variability of the random component and not the one of the "raw" data. Or am I missing something? Has anybody stumbled upon a similar problem before? How would you guys go about it? As always, any help is greatly appreciated!


title: CDF of largest observation in normal distribution

text: Let $X_1,...,X_n$ be a random sample from a $\mathcal{N}(\mu,1)$ distribution. Only the largest observation $Y = \max(X_1,...,X_n)$ is reported. What is the density of $Y$? How do I get there?


title: How to test for changes in 6 yes / no questions assessed before and after treatment

text: Young children were given a questionnaire before and after a presentation. There were 6 questions to answer yes or no. I would like to find out if there is a significant difference between the questionnaires taken before and after the presentation. I could score each test and do repeated measures testing, but the questions do not necessarily have right or wrong answers, so I would like to treat it as binary data. I would also like to test for each question. For example, on the first question, there were 12 yes answers before the presentation, and 43 yes answers after.


title: Why is multiple comparison a problem?

text: I find it hard to understand what really is the issue with _multiple comparisons_. With a simple analogy, it is said that a person who will make many decisions will make many mistakes. So very conservative precaution is applied, like Bonferroni correction, so as to make the probability that, this person will make any mistake at all, as low as possible. But why do we care about whether the person has made any mistake at all among all decisions he/she has made, rather than the _percentage_ of the wrong decisions? Let me try to explain what confuses me with another analogy. Suppose there are two judges, one is 60 years old, and the other is 20 years old. Then Bonferroni correction tells the one which is 20 years old to be as conservative as possible, in deciding for execution, because he will work for many more years as a judge, will make many more decisions, so he has to be careful. But the one at 60 years old will possibly retire soon, will make fewer decisions, so he can be more careless compared to the other. But actually, both judges should be equally careful or conservative, regardless of the total number of decisions they will make. I think this analogy more or less translates to the real problems where Bonferroni correction is applied, which I find counterintuitive.


title: How to identify transfer functions in a time series regression forecasting model?

text: I am trying to build a time series regression forecasting model for an outcome variable, in dollar amount, in terms of other predictors/input variables and autocorrelated errors. This kind of model is also called dynamic regression model. I need to learn how to identify transfer functions for each predictor and would love to hear from you about ways to do just that.


title: How to measure a classifier's performance when close to 100% of the class labels belong to one class?

text: In my data, I have a class variable, denoted as $C$. This class variable values are ${0, 1}$ (binary). Almost all observations of $C$ are 0 (close to 100%, more precisely, 97%). I would like a "performance" test on different classification models (it could be accuracy). What I am afraid of happening is that if I have a classification model that always classifies any observation into class 0, then that model will be 97% accurate (even though it never considers any other variables). Are there any well known performance tests for classification models on data dealing with very rare events?


Do you ever center AND standardize variables in multiple regression? It seems as if standardization would automatically center variables...is this true?


title: Performance of the fuzzy c-means clustering algorithm

text: I have implemented a genetic algorithm for a fuzzy c-means clustering in Matlab. Its performance should be apriori better than that of the classic fuzzy c-means (fcm function in matlab). However, on all the classical datasets I am comparing them, both algorithms converge to the same value of the objective function. My question is the following: is it just my bad luck, or FCM usually converges to the global minimum? Are there any datasets known when FCM shows suboptimal performance? Thanks,


title: Wavelet smoothing at different scales

text: I would like to perform wavelet smoothing at different scales in R. I got the idea from this figure (panel A) were they measure the density/intensity of a particular signal at different scales. Here the image shows the linear coordinates along one chromosome on the x axis and the y axis shows the 'bandwidth' or bin size in which the signal was measured. Does using wavelets here makes sense? The first figure comes from de Wit et al. Global Chromatin Domain Organization of the Drosophila Genome, PLoS Genetics 2008. ![Taken from "Global Chromatin Domain Organization of the Genome"](http://i.stack.imgur.com/3fqMo.jpg) I found this kind of multi-scale wavelet plot but I am not sure if it is reflecting the same thing as the first figure. This figures comes from this website ![This image shoes data from tree rings](http://i.stack.imgur.com/Ez9OO.jpg)


title: Estimating 2D mean from uneven samples

text: This CrossValidated post on estimating mean from uneven samples suggests using a (any) Riemann Sum to estimate the integral of an unevenly sampled 1D function. To me the obvious extension to 2D is this: given an uneven sample of f(x,y), create a triangulation of your sample points (x_i,y_i,z_i). Estimate the integral by adding the volume under each triangle. However, triangulating 3D points is not a trivial problem, so I imagine your results might depend on the quality of your triangulation. Is there a better/simpler way of creating a piecewise planar approximation to a 2D surface if not in general, at least for the purpose of computing the integral?


title: Which pareto distribution?

text: How do I assess what Pareto (or similar) function/distribution is most appropriate for fitting my data? I want to be able to model a dataset with a distribution of many-small, few- large that looks very similar to a Pareto distribution, but I do not have the statistics background to evaluate what model is most appropriate.


title: How to find the Likelihood function of superposed processes for MCMC?

text: i'm trying to understand how MCMC works. In my special case I have to stochastic processes which are superposed. I observe data $Y=\\{Y_1,...Y_T\\}$ and assume, that there are some latent variables $A_t$ and $B_t$ such that $Y_t = A_t + B_t$ for all $t$. Now, $A$ is generated by an Ornstein Uhlenbeck Process while $B$ is generated by a standard brownian motion. I managed to implement a MCMC simulation study using Random Metropolis Hastings for both processes individually, but I don't know how to combine them. The single likelihood function (or transition density) for the OU process is given by $p(t,x,y) = \frac{1}{\sqrt{\pi \sigma^2(1-\exp\left[ -2 \beta t\right])/\beta}}\exp \left[ - \frac{(y-\alpha-(x-\alpha) \exp{- \beta t})^2}{\sigma^2 ( 1 - \exp \left[ -2 \beta t \right])/\beta} \right]$ and the likelihood function of the brownian motion is given by $p(t,x,y) = \frac{1}{\sigma \sqrt{2 \pi t}} \exp \left[- \frac{1}{2 \sigma^2 t} (y-x-\mu t)^2\right]$. Is there an analytical way to get the joint likelihood function or do I have to estimate it somehow numerically? If so, some advice would be really appreciated. I also read something about the Clifford Hammersley theorem, which might be a solution for my problem. But I don't understand how to apply it in my special case. (How to draw from the conditionals?!) Thanks!


title: Inference on bootstrapped-confidence interval resulting in multimodal distribution

text: I have a set of paired data that I am using to calculate a ratio. The ratio is calculated as the sum of column 2 over the sum of column 1 (please see below)               Col1    Col2     3744519.37  0     2739505.37  0     46378122.62 0     2887639.57  78652.84     6000000 0     971112.48   0     7815969.01  0     11920577.36 0     2060120.66  0     4951863.34  0     8559269.86  0     4246529.62  0     54250000    0     17893683.49 0     1777305.52  0     1800832.37  0     2838034.29  0     15291266.83 0     3096538.88  0     2088767.28  0     5981199.57  0     2141878.34  0     1684058.06  0     2694327.27  0     7736276.87  0     3356067.88  0     1819461.59  0     972348.2    0     6922365 0     3357005.59  0     14870965.09 0     21500000    0     22419135    0     10000000    0     9250007.22  0     8900000 0     10404742.45 0     20452123.18 0     6340460.82  0     30494423.14 0     1157659.68  0     2868406.71  0     94257879.71 0     10369155.12 0     4315000 0     6158880 0     20057539.59 0     1335834.38  0     2406177.51  0     905846.54   405847      What I want to do is obtain a confidence interval from this point estimate to determine what the range of the population's ratio is. I took the approach of doing a bootstrap in R (see below)               iterations = 10000          PCarray = c()          for (i in 1:iterations){          rows=sample(1:length(data$Col1), length(data$Col1), replace = TRUE)          numerator   =  0     denominator =  0     ratio = 0          for (j in 1:length(rows)){          numerator   =  numerator   +  data$Col2[rows[j]]         denominator =  denominator +  data$Col1[rows[j]]          }          ratio = numerator / denominator          PCarray = c(PCarray,ratio)          }          hist(PCarray,breaks=200)      Which ultimately gives me this graph: ![Bootstrap](http://i.stack.imgur.com/f28Hl.png) As before, I want calculate a confidence interval. Is it alright just to do the percentile method to get a confidence interval (i.e. the upperbound is the 95th quantile and the lowerbound is the 5th quantile for a 90% CI)? I realise the driver in the results is that only 2 values are non-zero, but is this an acceptable approach? Is there any literature on this that states a confidence interval can be drawn from a non-symmetric, multimodal distribution?


title: Comparison of time series data between groups of individuals

text: We are currently analyzing experimental data with the following format: All participants (n=100) make repeated allocation decisions in fixed intervals: 100% is divided up into two allocation options at each step, therefore the dependent variable can be treated as a proportion. For each participant we therefore obtain a time series of proportions. In addition, there are two between-subjects factors with two levels each. We would like to find an appropriate model, in which can test the main effects and interactions of the two factors (or ideally three factors, taking time points into account). It would be nice to integrate components of ARIMA-models (the information obtained at different points in time make it likely that all participants will show systematic fluctuations regardless of group). Does anyone knwo appropriate models thart could be used here (groups of time series in a factorial design with proportions as dependent variable)? There is an additional complication (we can restrict ourselves to the problem sketched above, but it would be interesting if we could integrate this additional part): all groups overlap for about 70 periods, one group starts already 100 periods earlier, a second grouop about 70 periods earlier. In other words: group start at periods 1, 30, 100 and 100 and we collect 170, 140, 70 and 70 observations respectively. Would any model be flexible enough to integrate the additional information for the first two groups for building the time series model and still be able to allow for a valid hypothesis test for the two experiemntal factors for the 70 shared observations only?


title: Which test is appropriate for comparing Pre and Post intervention data?

text: I am conducting research in which employee engagement will be evaluated on approximately 20 participants. They will complete an 11 question Likert-scale survey, have interventions done on them (read a book, read 3 articles, and sit through a lecture), and complete the same 11 question survey. Is the correlated t-test the appropriate test to conduct on the pre and post data to determine if a statistically significant increase in employee engagement occurred?


title: Alternative for mi impute intreg

text: Is there a free alternative available for the Stata procedure mi impute intreg (Impute using interval regression)? For example as an R package. I have not found any yet.


title: A multi-label classification for tagging short text

text: I am fairly new in the area of text mining and want to practice my skills a little. I have the following task at hand which I want to work on. I have a large list of short texts (~100.000) and every text has on average 3 tags / labels assigned to it. What I want now is train a predictive model, with which I can assign labels to new unseen texts. The assumption about the labels is, all existing labels are included in my dataset. I started with python sklearn and a basic tf.idf matrix representation combined with an one-vs.-rest support vector classification for each label, but the mean f1 score for that model are very low. Are there better ways of creating a predictive model for that kind of task? As I mentioned before, this area is quite new to me so a little advice on how to tackle this problem would be more than welcome (e.g. some literature ?).


title: How to perform a trend test

text: I'am working on a glm model (family:normal, link:identity) using complex survey design (package: survey), My interest is to evaluate the effect of the variable DLA on the scores. I ran the model 1 considering DLA as categorical (tertiles):               model1=svyglm(score~DLA_tert+race+pir_4+age_yr+as.factor(educa),design1[age_yr>69] ,data=base_sex2[age_yr>69])      Here are the results for this model               Coefficients:                         Estimate Std. Error t value Pr(>|t|)         (Intercept)          90.6470    13.4223   6.753 6.49e-06 ***     DLA_tert(37.1,63.1]  -1.5940     1.9428  -0.820 0.424793         DLA_tert(63.1,507]   -3.1713     2.2090  -1.436 0.171632         race2               -12.6004     2.4078  -5.233 0.000101 ***     race3                -9.1128     2.5326  -3.598 0.002635 **      race4               -10.3232     2.4077  -4.288 0.000648 ***     pir_42                2.7497     2.5701   1.070 0.301593         pir_43                6.8206     2.6690   2.555 0.021959 *       pir_44                6.3329     3.2067   1.975 0.066985 .       age_yr               -0.8075     0.1779  -4.539 0.000392 ***     as.factor(educa)2     8.3018     2.3218   3.576 0.002760 **      as.factor(educa)3    13.5644     2.0637   6.573 8.84e-06 ***     as.factor(educa)4    15.3992     2.2397   6.876 5.27e-06 ***     as.factor(educa)5    17.9402     2.7926   6.424 1.14e-05 ***          I suspect there is a trend in the association between my tertiles and scores (even if the association is not significant). I would like to know how to calculate this trend and the p-value for it? Thank you very much


title: Regularized bayesian logistic regression in JAGS

text: There are several math-heavy papers that describe the Bayesian Lasso, but I want tested, correct JAGS code that I can use. Could someone post sample BUGS / JAGS code that implements regularized logistic regression? Any scheme (L1, L2, Elasticnet) would be great, but Lasso is preferred. I also wonder if there are interesting alternative implementation strategies.


title: Real utility of small accuracy improvements in sentiment classifiers

text: I have lately been reading papers regarding Sentiment Analysis, where most researches report that their improvements made them achieve an increase of 1~2%, or even 0.5% in accuracy compared to baseline (non-trivial) methods. Of course, I understand that such an increase (even if small, if statistically significant) is a good thing, but what are the advantages, in terms of application/utility of such a small improvement ? Is there any application/practical use that could directly benefit from such a small improvement ? If such an application exists, how to know which is the minimal accuracy needed for it? (Similar question is posted here, but I am asking for a more practical view rather than statistical) Thanks!


title: How to run two-way ANOVA on data with neither normality nor equality of variance in R?

text: I am working on my master thesis at the moment and planned on running the statistics with SigmaPlot. However, after spending some time with my data I came to the conclusion that SigmaPlot might not be fit for my problem (I may be mistaken) so I started my first attempts in R, which did not exactly make it easier. The plan was to run a simple TWO-WAY-ANOVA on my data which results from 3 different proteins and 8 different treatments on those, so my two factors are proteins and treatments. I tested for normality using both               > shapiro.test(time)      and               > ks.test(time, "norm", mean=mean(time), sd=sqrt(var(time)))      In both cases (maybe not surprising) I ended up with a non-normal distribution. Which left me with the first questions of which test to use for equality of variances. I came up with               > chisq.test(time)      and the result was, that I don't have equality of variance in my data either. I tried different data transformations (log, center, standardization), all of which did not solve my problems with the variances. Now I am at a loss, how to conduct the ANOVA for testing which proteins and which treatments differ significantly from each other. I found something about a Kruskal-Walis-Test, but only for one factor (?). I also found things about ranking or randamization, but not yet how to implement those techniques in R. Does anyone have a suggestion what I should do? Edit: thank you for your answers, I am a little overwhelmed by the reading (it just seems getting more and more instead of less), but I will of course keep going. Here an example of my data, as suggested (I am very sorry for the format, I couldn't figure out another solution or place to put a file. I am still new to this all.):               protein treatment   time       A   con 2329.0       A   HY  1072.0       A   CL1 4435.0       A   CL2 2971.0       A   CL1-HY sim  823.5       A   CL2-HY sim  491.5       A   CL1+HY mix  2510.5       A   CL2+HY mix  2484.5       A   con 2454.0       A   HY  1180.5       A   CL1 3249.7       A   CL2 2106.7       A   CL1-HY sim  993.0       A   CL2-HY sim  817.5       A   CL1+HY mix  1981.0       A   CL2+HY mix  2687.5       B   con 1482.0       B   HY  2084.7       B   CL1 1498.0       B   CL2 1258.5       B   CL1-HY sim  1795.7       B   CL2-HY sim  1804.5       B   CL1+HY mix  1633.0       B   CL2+HY mix  1416.3       B   con 1339.0       B   HY  2119.0       B   CL1 1093.3       B   CL2 1026.5       B   CL1-HY sim  2315.5       B   CL2-HY sim  2048.5       B   CL1+HY mix  1465.0       B   CL2+HY mix  2334.5       C   con 1614.8       C   HY  1525.5       C   CL1 426.3       C   CL2 1192.0       C   CL1-HY sim  1546.0       C   CL2-HY sim  874.5       C   CL1+HY mix  1386.0       C   CL2+HY mix  364.5       C   con 1907.5       C   HY  1152.5       C   CL1 639.7       C   CL2 1306.5       C   CL1-HY sim  1515.0       C   CL2-HY sim  1251.0       C   CL1+HY mix  1350.5       C   CL2+HY mix  1230.5


title: Maximizing: likelihood vs likelihood ratio

text: Say I have an **observed** data set ($n_i$) and I want to obtain the best fit out of 10 data sets produced by a **model** dependent on a single parameter $a$ ($m_i(a)\;a=1..10$). Suppose I use a Poisson likelihood distribution: $P_i(a)=\frac{m_i(a)^{n_i}}{e^{m_i(a)}n_i!}\; ; \; a=1..10$ where $m_i(a)$ represents the **model** value of bin $i$ (for a given value of the parameter $a$) and $n_i$ its **observed** value (which remains constant always) The likelihood of the whole data set or _cumulative likelihood_ (for each value of $a$) is then: $L(a)=\prod\limits^{N}_{i=1} \frac{m_i(a)^{n_i}}{e^{m_i(a)}n_i!}\; ; \; a=1..10$ where $N$ is the total number of bins. Now, I want to pick the **model data set** that best fits my **observed data set**. Could I just choose the maximum value of $L(a)$, ie: the _maximum likelihood_ : $Max\\_Likelihood = max\,[ L(a)\; ; \; a=1..10]$ and say that the value of parameter $a$ associated with that particular **modeled data set** is the best estimation of parameter $a$, given my **observed data set**? Or should I use a _cumulative likelihood **ratio**_ , defined as : $LR(a)= \prod\limits^{N}_{i=1} \frac{\frac{m_i(a)^{n_i}}{e^{m_i(a)}n_i!}}{\frac{n_i^{n_i}}{e^{n_i}n_i!}} = \prod\limits^{N}_{i=1} \left(\frac{m_i(a)}{n_i}\right)^{n_i}e^{n_i-m_i(a)}\; ; \; a=1..10$ and keep the value of $a$ that gives the maximum value of $LR(a)$, ie the _maximum likelihood **ratio**_ : $Max\\_Likelihood\\_Ratio = max\,[ LR(a)\; ; \; a=1..10]$ Since the observed data set ($n_i$) remains the same for all the modeled data sets ($m_i(a);\;a=1..10$), won't maximizing this value $LR(a)$ give me the same result (ie: the same value of $a$) as maximizing $L(a)$?


title: Optim error - initial value in vmmin is not finite

text: I have run                 'fitdistr(mydata, "Gaussian")'       and get the following error message:               Error in optim(x = c(34.513, 25.587, 18.251, 52.836, 38.206, 42.395, 47.22,  :        initial value in 'vmmin' is not finite      However when I run                'fitdistr(mydata,"Normal")'      I get the following:                 mean           sd           39.71581691   23.44081474      ( 0.02856006) ( 0.02019501)      How do I fit my data to Gamma and other distributions that are not normal? Sajni


title: Joint confidence intervals for probabilities

text: I have two probabilities $p$ and $q$. $p>q$, and they aren't correlated. I'm going to calculate $i$ such that $p^i=q$, which is easily done as $\log_p(q)$. Now, I'd like to also calculate a confidence interval for $i$, which is necessarily going to be a function of both $p$ and $q$'s confidence intervals. My first approach was to do               p.min <- qbeta(0.025, 152, 29)     p.max <- qbeta(0.975, 152, 29)          q.min <- qbeta(0.025, 37, 19)     q.max <- qbeta(0.975, 37, 19)          ## q.max and p.min have the smallest difference     i.min <- log(q.max, base = p.min)          ## q.min and p.max have the largest difference     i.max <- log(q.min, base = p.max)      But it occurs to me that 95% confidence intervals for $p$ and $q$ independently probably produces too large a confidence interval for $i$, because the joint confidence interval for $p$ and $q$ will be narrower. So, how do I go about figuring out the joint confidence interval of $p$ and $q$. They're uncorrelated, which should make things easier. It is as simple as narrowing the quantiles in `qbeta()`? By how much?


title: Why one can fit Poisson GLM on non-integer values in Stata?

text: I have aggregated death counts and some categorical predictors e.g. era, drug use, and I'm running a Poisson GZLM in SPSS, also a `glm` in Stata (generalized linear model). There are many, but not excessive, zeros and it has been suggested to me to add 0.5 to the observed values. When I do that the GENLIN procedure in SPSS won't run, because the Poisson dependent should have integer values. However, the `glm` with family Poisson and log link in Stata runs and just makes a note that the dependent has non-integer values. Why does Stata run the Poisson `glm` if the values of the dependent should be integers only?


title: What is the standard format of the equation that represents the production of this daughter isotope?

text: I have an idea about producing a genetic algorithm. But it involves using a relationship which believe could be defined as a horizontal Asymptote. I have found an image that depicts the type of relationship I will attempt to apply. However finding out what it is called and what the standard equation is that defines it is escaping me. # Decay relationship of interest ![a not so busy cat](http://i.imgur.com/wdQgUBW.gif) I am interested in what the standard format for the green dashed graph would be. It was a long time since I took any calculus. So any help would be appreciated. AEA


title: How to use liblbfgs for fitting?

text: I am trying to use bfgs algorithm in order to fit a set of $\\{(x,y),f(x,y)\\}$ to a function in the form of let's say $a\cdot cos(x)+b \cdot y=f(x,y)$. I try to understand how to use bfgs algorithm with liblbfgs, but I don't understand the example and it is not clear what function the author tried to fit.


title: Fitting a gam model with simple Gaussian/Student-t heteroskedasticity

text: I am fitting a gam model in R (using the `gam` function in `mgcv`) to account for some non-linear effects in my data. A stripped down example of what I am doing in R is:               mod=gam(y~s(x)+s(z),data=df)      However, I want to add a slightly more complicated variance model to my regression of the form $$\epsilon \sim N(0,\sigma^2),\ \sigma = f(\hat{\mu})$$ where $\hat{\mu}$ is the fitted value of the model. (Actually, it would be nice if $\epsilon \sim t_\nu$ for some $\nu$ but sticking to this for now. I have managed to do this in the `gls` function from `nlme` using the `varFunc(form=fitted(.))` type approach, but can't figure out if there is an option to do the same kind of thing using `gam`. I recognise this is not really the intention of a GLM/GAM model, but I don't want to reinvent the wheel if I am just missing something obvious Edit: In response to the question in the comment below, I am hoping to fit a linear or quadratic function for $f$. I do not know the exact form of $f$ but plan to iteratively estimate it from the residuals if this can't be done automatically. Edit2: Typo in R code - first spline is not meant to be a function of y!


title: Persistent Contrastive Divergence for RBMs

text: When using the persistent CD learning algorithm for Restricted Bolzmann Machines, we start our Gibbs sampling chain in the first iteration at a data point, but contrary to normal CD, in following iterations we don't start over our chain. Instead we start where the Gibbs sampling chain in the previous iteration ended. In the normal CD algorithm each iteration evaluates a mini batch of data points and computes the Gibbs sampling chains starting from those data points themselves. In persistent CD, should we keep Gibbs sampling chains for _each_ data point? Or should we keep also a mini batch of Gibbs sampling chains, which started at data points which aren't currently evaluated in the current iteration? It seems to me that keeping Gibbs sampling chains for each data point will be too cumbersome, but on the other hand it seems inadequate to compare the signals of the current sample with the signals after a long Gibbs chain _which didn't start at the current sample_.


title: Denominator is Zero for Matthews correlation coefficient and F-measure

text: Recently, I built a classification model based on the imbalanced data set(positive sample is minority and negative sample is majority), and the model gave the following result for the test set: > True Positives = 0 > > True Negatives = 139 > > False Positives = 0 > > False Negatives = 10. My question is: for the result, can Matthews correlation coefficient (MCC ) and F-measure be used for estimating the classifier? Since the denominators for MCC and F-measure are zero, it seems meaningless. If so, MCC and F-measure is not always works for estimating the classifier, and sensitivity and specificity as well as g-mean should be better. Is that right? Any help is appreciated.


title: Residual structure in Linear Mixed Models with Random effects

text: You can further improve a linear mixed model with random intercept and slope by specifying a structure in the residuals (for example AR(1)). In SAS it is possible, but I hope this is also already possible in R. Does somebody knows how to do this? I am using both packages lme4 and nlme. Thank you, Kasper


title: Forecasting an ADL model

text: I am fairly new to Stata, currently taking an undergrad time series econometric class. The economic significance of the regression I am running and attempting to forecast is all but zero; this is simply my learning how to use Stata somewhat efficiently. I have run an Autoregressive Distributed Lag (ADL) regression of Industrial Production (`ip`) using two lags of `ip` and one lag of the Federal Funds Rate (`ffr`). My goal: forecast this model 30 periods into the future. I was hoping I could run an AR(2) model of the `ffr` and forecast it 30 periods (this is straightforward and I have done so) and then use these values to forecast my ADL model 30 periods. Does this make sense and could someone inform me on how I could do this or suggest another method?


title: Aligned Rank Transform with random factor

text: In an aligned rank transform preceding a two-way ANOVA, the procedure is:   * save residuals by performing a standard ANOVA    * use Aggregate to determine effects for group means (mij for interaction, ai as first factor, bj for second factor)    * eliminate interaction effect from residual to determine interaction (ran), eliminate single effects to determine main effects (ra)  Interaction effect: ek – (αβij – αi – βj + μ) Main effects: ek – (αi + βj–2μ)   * transform residuals into ranks    * run complete model on these ranks for interaction / main effects, respectively My question: how to incorporate a third RANDOM factor in addition to these two fixed factors ai and bj? Will it be included just as the others, or will random factors have to be treated differently? Thanks!


title: Using chi goodness of fit to test if data comes from normal distribution

text: The problem gives 50 observations and then asks to check whether the data comes from a normal distribution. For instance if these are your observations then how would you go about testing to see if they come from a normal distribution.               −1.28 −1.22 −0.45 −0.35 0.72     −0.32 −0.80 −1.66 1.39 0.38     −1.38 −1.26 0.49 −0.14 −0.85     2.33 −0.34 −1.96 −0.64 −1.32     −1.14 0.64 3.44 −1.67 0.85     0.41 −0.01 0.67 −1.13 −0.41     −0.49 0.36 −1.24 −0.04 −0.11     1.05 0.04 0.76 0.61 −2.04     0.35 2.82 −0.46 −0.63 −1.61     0.64 0.56 −0.11 0.13 −1.81      The professor talked something about transforming it into a uniform distribution and then having 4 intervals 0,0.25, 0.5,0.75,1 and then see how many values fall within those intervals but I'm not sure how to go about doing that.


title: Simple approximation of Poisson cumulative distribution in long tail?

text: I want to decide the capacity $C$ of a table so that it has residual odds less than $2^{-p}$ to overflow for given $p\in[40\dots 120]$, assuming the number of entries follows a Poisson law with a given expectancy $E\in[10^3\dots 10^{12}]$. Ideally, I want the lowest integer `C` such that `1-CDF[PoissonDistribution[E],C] < 2^-p` for given `p` and `E`; but I'm content with some `C` slightly higher than that. Mathematica is fine for manual computation, but I would like to compute `C` from `p` and `E` at compile time, which limits me to 64-bit integer arithmetic. Update: In Mathematica (version 7) `e = 1000; p = 40; c = Quantile[PoissonDistribution[e], 1 - 2^-p]` is `1231` and seems about right (thanks @Procrastinator); however the result for both `p = 50` and `p = 60` is `1250`, which is wrong on the unsafe side (and matters: my experiment repeats like $2^{25}$ times or more, and I want demonstrably less than $2^{-30}$ overall odds of failure). I want some crude but safe approximation **using 64-bit integer arithmetic only** , as available in C(++) at compile time.


title: GAM sample size

text: A common rule of thumb in social science is to have at least $30$ observations for regression (without any complicated effects). Does someone know such rule of thumb for GAMs? Are $50$ observations enough for a simple additive model with two parametric terms and two smoothers?


title: Cross validating logistic GAM using CVgam

text: I have an ecological data set, whereby sediment mud content (%) (i.e. the continuous explanatory variable) is thought to be explaining the spatial distribution (i.e. presence/absence) of various estuarine benthic macroinvertebrate taxa (i.e. the binary response variable). Below is a subset of the data (actual data consists of >1400 replicates) for a single taxon, where the response is 'Bin' and the single predictor is 'Mud' and 'Estuary' is considered a random effect.               Estuary       Mud   Bin     FW            0.1   0     FW            0.3   0     FW            0.2   0     FW            2.1   0     POR           7.1   0     FOR           11.4  0     JRE           25.5  1     JRE           33.6  1     JRE           44.6  1     JRE           55.6  1     JRE           63.6  1     FW            76.6  1     FW            87.6  1     FW            90.6  1     FW            93.6  1      Here's the GAM model which used the `mgcv` package (and was conditioned on the 1400+ data points):               aa1.estuary<-gam(Bin~s(Mud,bs="ps", k=6) + s(Estuary, bs="re"),family=binomial, gamma=1,data=sample)      Here's the Cross validation step, which used `CVgam` in the `gamclass` package, but returned the error below when `Estuary` was included as a random effect:               aa1.val <- CVgam(formula=Bin~s(Mud, bs="ps", k=6) + s(Estuary, bs="re"), data = sample, nfold = 10, debug.level = 0, printit = TRUE, method = "GCV.Cp",cvparts = NULL, gamma = 1, seed =100)     Error in X[ind, ] : subscript out of bounds      To check, `Estuary` was omitted, generating the following output:               aa1.val <- CVgam(formula=Bin~s(Mud, bs="ps", k=6) , data = sample, nfold = 10, debug.level = 0, printit = TRUE, method = "GCV.Cp",cvparts = NULL, gamma = 1, seed =100)          GAMscale    CV-mse-GAM            0.1491      0.1499      Hence I have two related questions: (1) whats causing the error when `Estuary` is included in the cross-validation step? and (2) do the above estimates of `GAMscale=0.1491` and `CV-mse-GAM =0.1499`seem too close to one another to be realistic? Thanks in advance,


title: Terminology and handling of log-normal mixture distributions

text: The definition of a _log-normal distribution_ of a random variable is based on normality of its logarithm. I'm curious whether there exist a _specific term_ for cases, where _log-transformed data_ does not conform to a normal distribution, but to a _mixture distribution_. Moreover, I would like to know what is the appropriate **procedure(s)** (both _visual_ and _analytic_ ) of _demonstrating_ the presence of a mixture distribution in log-transformed data by superimposing a plot of the log-transformed data with its component distributions (probably, using _kernel density_ plots). With that in mind, I'm interested in applied solution in `R`, preferably using `ggplot2` package. This question is based on my attempts to implement the above-mentioned demonstration for my _exploratory data analysis_. For a _normal mixture distribution_ , I have figured out the solution (http://stackoverflow.com/a/25588144/2872891), but I'm still struggling with _log-normal mixture distribution_ (http://stackoverflow.com/q/25598485/2872891). Any help or advice is welcome!


title: How do you identify the class of an anonymized attribute in statistical analysis?

text: How do you identify the class of an "anonymized attribute" in statistical analysis? (By "anonymized attribute" I mean that I do not know whether the attribute is continuous or categorical). The attribute takes integer value. I would like to know whether the attribute is continuous or categorical (nominal or ordinal). How I can identify that?


title: How to calculate the prediction interval from two independent forecasts of time series?

text: Lets say I have two time series a and b and I predict the next value and a prediction interval (say lower and upper 80%) of each of these two. I end up with something like that: a: [10, 20, 30] and b: [2, 7, 12] (as [lo80%, prediction, high80%]). As these time series are independent, I see I could just add 20 + 7 to get a prediction of a + b, but what about the prediction interval? I suppose I cant just add them to [12, 27, 42], right? (Sorry if it is too simple a question)


title: GARMA models for counts

text: Suppose we have count data. It seems that a generalized ARMA model would be good to use. In Rm these is function: > garma(link = c("identity", "loge", "reciprocal", "logit", "probit", > "cloglog", "cauchit"), earg=list(), p.ar.lag = 1, q.lag.ma = 0, coefstart = > NULL, step = 1) Why type of link function would be good to use? There is no poisson option.


title: Can I use a one way ANOVA with multiple categoricals as a response variable?

text: I do not know a lot about stats so forgive me. I have a dataset which has $N$ number of participants and their choice between 1, 2, and 3 at any given point. I want to analyze the distribution of the choices between 1, 2, and 3 by participants' age and sex. How would I set up a test to look for significant interactions between these variables? I thought I needed a one way ANOVA but I am not sure. Thank you


title: How to perform k-means clustering with only a distance function, not euclidean points?

text: I want to perform k-means clustering on some objects I have, but the objects aren't described by "points". However, I am able to compute the distance between any two objects (it is based on a similarity function). I've implemented K-means before, but it's not clear to me how to update the clusters to be the cluster "centers" without a point-representation. How would this normally be done?


title: Filtering techniques and noise

text: Suppose we have some house price data for 30 years (1970-1999). This is yearly data (30 data points). Suppose some major event $X$ happened on 1980. I want to see whether this event affected prices later on. Is using median filters a good methodology for detecting outliers in later years?


title: How are ranks calculated in tie cases?

text: I am considering the following question: when there are ties, how are ranks calculated? In some references, they first rank them without repeating the ranks, and then average the ranks of those ties and assign the average rank to each one in the ties. I was wondering if that way is unanimously used in statistics for determining the ranks? Are there other ways to determine the ranks in tie cases? If yes, when to use which way? Thanks and regards!


title: Cosine similarity for a group of vectors

text: I am starting to work in R with the following case: 35 countries have scores between [0,1] in 120 dimensions and I need to show that there are differences in the way their scores distribute across these dimensions. For example, if there were 5 dimensions, (1,0.5,0,0,0) is different than (0,0,0,0.5,1). This is why I chose to use Cosine Similarity (CS) as a measure, since I am more worried about the "direction" rather than the intensity of their scores. Since CS is a measure between individual vectors, I end up with a 35x35 matrix of values but my problem is that I still need to find an appropriate criteria/test to say if they are **all** similar or not. I'm uncertain if it would be correct to just say that if all values in the matrix are >0.95 then they are all similar. I'm also not very good with statistics, so forgive me if I said some nonsense. Any advice? PD: I also thought about seeing if data is actually clustered or not, i'm not sure.


title: Monte Carlo Rejection Sampling Method

text: I have the following passage from a set of lecture notes I am working on that I would like to understand a little better. $\underline{\text{Algorithm for Rejection Sampling}}$: Given two densities $f$, $g$, with $f(x) < M.g(x)$ for all $x$, and some constant $M$, we can generate a sample from $f$ by   1. Draw $X \sim g$   2. Accept $X$ as a sample from $f$ with probability $\large \frac{f(X)}{M.g(X)}$ otherwise go back to step 1. Now, I understand why this will generate samples from the distribution with density $f$, but I really do not understand how to use this algorithm in practice. I draw an $X$ from the distribution with density $g$, and then I need to decide whether to keep it, or reject it. I don't understand the condition to accept it... do I accept it if the probability defined in 2. is greater than e.g. 0.5, or maybe 0.7? Is this probability up to me to decide? Thanks for your insight and thoughts.


title: Is Wikipedia wrong about the null distribution of the signed rank statistic?

text: From Wikipedia > Calculate the test statistic $W$ $$ W = \left|\sum_{i=1}^{N_r} > [\operatorname{sgn}(x_{2,i} - x_{1,i}) \cdot R_i]\right|, $$the absolute > value of the sum of the signed ranks. > > As $N_r$ increases, the sampling distribution of $W$ converges to a normal > distribution.   1. Since $W$ is defined to be always nonnegative, why does the distribution of $W$ converge to a normal distribution which has positive probability of negative values? What should it be then?   2. > For $N_r \ge 10$, a z-score can be calculated as $z = \frac{W - 0.5}{\sigma_W}, \sigma_W = \sqrt{\frac{N_r(N_r + 1)(2N_r + 1)}{6}}$. Does it mean $E(W) = 0.5$ and $\operatorname{Var}(W) = \frac{N_r(N_r + 1)(2N_r + 1)}{6}$? Why is it true? Thanks and regards!


title: Getting Residuals to be White Noise

text: I'm on a time series project for an undergraduate course. For the project I'm trying to come up with an ARIMA model for the housing starts data set. http://www.quandl.com/FRED/HOUST-Housing-Starts-Total-New-Privately-Owned- Housing-Units-Started I have been trying to adjust the p and q to get residuals that are uncorrelated white noise. I tried auto.arima to get parameter values.               > fit.start=auto.arima(data)     > fit.start     Series: data      ARIMA(2,1,3)(2,0,0)[12]                              Coefficients:              ar1      ar2      ma1     ma2      ma3     sar1     sar2           1.2874  -0.3256  -1.6559  0.8715  -0.1725  -0.1352  -0.2074     s.e.  0.5709   0.5348   0.5665  0.7266   0.2052   0.0520   0.0517               > fit213200<-Arima(data, order=c(2,1,3), seasonal=c(2,0,0))     > tsdisplay(residuals(fit213200))      ![enter image description here](http://i.stack.imgur.com/sBbjL.png) So is there anyway to get the residuals to be just white noise?


title: Machine Learning Algorithms vs. Linear Regression

text: Do machine learning algorithms like Boosted Regression Trees (in the R package (gbm)) follow the same statistical assumptions of not including correlated predictor variables in GLM? i.e. If I have two correlated predictrs (rsq=.7) should I be including both into my BRT model? Any input or thoughts on this question would be greatly appreciated. Cheers,


title: Name for a function mapping time slices to probabilities

text: Is there a commonly used name for a function that maps time slices in a day to probability of a specific event happening in that slice? I was using probability distribution but I guess since the total of y values can be as large as the number of time slices (i.e., this event happens at every time slice) this is misusing the terminology.


title: Document classification question

text: I've got a tricky document classification question. Each document in my collection can be assigned a label. In my training set the average is about 3 different labels assigned to each document. What I have learned in the training set is $P(\text{label}|\text{doc_word})$ for many label/doc_word pairs. I want to calculate $P(\text{label}|\text{word}_1,\ldots,\text{word}_n)$ but I cannot simply calculate it by $$P(\text{label}|\text{word}_1,\ldots,\text{word}_n)=P(\text{label}|\text{word}_1)\cdot \ldots \cdot P(\text{label}|\text{word}_n)$$ because there are many words where $P(\text{label}|\text{word}_i)=0$. Furthermore, if there are two words in the document where $P(\text{label}|\text{word})>0$, their product would be less than the individual probabilities thus making the label less likely. This is counter- intuitive to me because I would think that having two words where $P(\text{label}|\text{word})>0$ should make the label more likely. So maybe the independence assumption does not hold. Perhaps I should just use these probabilities as association rules in order to assign labels?


title: Obtaining marginal distributions from the bivariate normal

text: Let $(X, Y)$ have a normal distribution with mean $(\mu_X, \mu_Y)$, variance $(\sigma_X^2, \sigma_Y^2)$ and correlation $\rho$. I want to know the corresponding marginal densities. All I found so far was the well-known density expressions for $X\sim N(\mu_X, \sigma_X^2)$ and $Y\sim N(\mu_Y, \sigma_Y^2)$, but isn't that just for $X \perp Y$? Shouldn't $\rho$ appear in the expressions? How do I calculate the marginals of any $X$ and $Y$? I think using the definition will end up in an integral that cannot be solved analytically…


title: Distributions of $X+Y$,$X-Y$,$XY$ for $(X,Y)$ Chosen Uniformly Inside Triangle

text: Let $(X,Y)$ be chosen uniformly on the triangle $\\{(x,y)\in\mathbb R^2:x+y\leq1,x\geq0,y\geq0\\}$. What is the density function of $(X,Y)$? Find the distributions of $X+Y$, $X-Y$,$XY$. What I've tried: $\displaystyle \frac{1}{Area\hspace{1mm} of\hspace{1mm} Triangle} =\frac{1}{\frac{1}{2}\cdot 1 \cdot 1}=2$, Therefore the joint density function of $(X,Y)$ is: $$ f_{X,Y}(x,y)= \begin{cases} 2 & \text{inside the triangle} \\\ 0 & \text{elsewhere} \end{cases} $$ How do I then find the distributions? Thank you!


title: How do you report significant contrast results for a MANOVA in APA style?

text: Do you just report that MANOVA significance test is significant, or do you have to include contrast estimates (SPSS `K MATRIX`)? In the latter case, what are recommended formatting guidelines?


title: Non-linear least squares SUR

text: I need to do a seemingly unrelated regression of a non-linear model. Can't figure out how to do it in Stata. Any advice much appreciated. Thanks, Ben


title: Linear regression from data that don't represent a function

text: I have $(x,\ y)$ pairs with a strongly suspected linear correlation. So I want to fit the "best" linear function in order to make predictions for unknown $x$'s. These pairs don't represent a function, i.e. there are many different values of $y$ for some $x$. I want to know if there is an issue with that considering this as a machine learning problem. Should I use some unique "representative" (average, maximum, minimum, most frequent, etc.) and convert the input data in a function or it is OK to work with the data as it is?


title: Using PC scores or cluster analsis in predictions

text: I have very big data and low number of observations. So I decided to use PCA to reduce dimension of the data. The following is R example (just an dummy example - for workout):               xmat <- matrix(sample(-1:1, 100000, replace = TRUE), ncol = 1000)     colnames(xmat) <- paste ("V", 1:1000, sep ="")     rownames(xmat) <- paste("S", 1:100, sep = "")      In this example dataset I have `1000` variables and `100` observations / subjects. I am doing PCA. Lets say.               out <- princomp(xmat)     Error in princomp.default(xmat) :        'princomp' can only be used with more units than variables      Q1: is there a way to reduce dimensionality with `p > n` ? I would like to use all variables information as opposed to representative ones. Without having proper solution I went anyway to use cluster analysis of variables to categorize the variables and pick the randomly from the clusters. To create a list of representative variables I tried to cluster the variables.               # cluster variables      d <- dist(t(xmat), method = "euclidean") # distance matrix     fit <- hclust(d, method="ward")     plot(fit)     groups = cutree(fit,40)     groupd <- data.frame(var = names(groups), group = groups)      What I am thinking is randomly pick one variable from each group above and use this in PCA. Assume that I have the following y variable.               set.seed(1234)     yvar.d <- data.frame (subject = c(paste("S", 1:100, sep = "")), yvar = rnorm (100, 50,10))      **Here is my question** :   1. What could be statistical challenge of using cluster analysis ?   2. Can we use PCA scores in predictions of y. How ? Just multiple regression or we can introduce something such as variance explained by each components in the model ? **Edits:** Based on the discussions (see the comments below), I am using different function to do PC analysis. "The calculation is done by a singular value decomposition of the (centered and possibly scaled) data matrix, not by using eigen on the covariance matrix. This is generally the preferred method for numerical accuracy. The print method for these objects prints the results in a nice format and the plot method produces a scree plot." - from function help.                    out1 <- prcomp(xmat)           out1$x[1:3,1:3]                           PC1        PC2       PC3     S1  2.940862 -2.7379835  6.527103     S2 -1.081124 -0.5294796 -0.276591     S3  2.375710  0.4505205 -4.236289             out1$sdev      screeplot(out1,npcs=30, type="lines",col=3) # 30 PCA plotted      ![enter image description here](http://i.stack.imgur.com/gMJys.jpg)                out1$rotation      I also come to see an example in SO how to use PCA in prediction. Here is my workout:               ## take our training and test sets     YY <-  yvar.d$yvar      prop <- 0.5     train = sample(1:length(YY), round(length(YY)*prop,0))               # data for testing model purpose      testid = setdiff (1:length(YY), train)     YY1 <- YY     newXPCA <- data.frame(out1$x)     test.data <- data.frame (y = YY1[testid],newXPCA[testid,])      test.data[1:10,1:10]          train.data <- data.frame(y= YY1[train],newXPCA [train,])     train.data[1:10,1:10]          ## fit the PCA     pc <- prcomp(train.data[, -1])     trainwPC <- data.frame (y = train.data$y, pc$x)          model1 <- lm(y ~ ., data = trainwPC)          #predict() method for class "prcomp"     test.p <- predict(pc, newdata = test.data)     pred <- predict(model1, newdata = data.frame(test.p), type = "response")     pred      Warning message:     In predict.lm(model1, newdata = data.frame(test.p), type = "response") :       prediction from a rank-deficient fit may be misleading      I just adopted this script from the SO link, I am not sure about accuracy of the script. I still have technical questions remaining such as clarification to **remaining question 2** above: (1) If I want to split data into training and test set by sampling `50% of data` (as show in the script). Should I do just multiple regression with y and the `out1$x` ? how many components to use ? is variance of each component play role in good model selection such as avoid over-fitting ? How ? (2) Clustering (using x clusters) vs PCA analysis (with subset of x components vs all ) what would be statistically favorite for predictions in the situations where have `p > n` ? As I said to my mind the PCA analysis can use all information but I do not know if there is downside of such information such as `over-fitting` and "error consumption". Worked example appreciated.


title: square of normal distribution with specific variance

text: What is the distribution of a normally distributed random variable $X^2$ with $X\sim N(0,\sigma^2/4)$. I know $\chi^2(1)=Z^2$ is a valid argument for when squaring a **standard normal** distribution, but in the case of non-unity variance?


title: Evaluating relative importance for a group of independent variables

text: I have done a multiple regression analysis of a group of design variables for a system I am analyzing. From this I determined the standardized regression coefficients. All of the variables fall into one of two groups, architectural variables and control variables. I would like to determine the relative importance of each group of variables, so that I can make a claim to the effect of "the architectural variables are twice as important as the control variables for this system". Is it acceptable if I simply sum the standardized regression coefficients from all the variables in each group and compare these two numbers?


title: Total Sum of Squares (TSS) - Exponential Regression

text: The total sum of squares is expressed as below [1][2], \begin{align} \rm{TSS} &= MSS + ESS \end{align} \begin{align} \sum_{i=1}^n (y_{i}-\overline{y})^2 &= \sum_{i=1}^n (\hat{y}_{i}-\overline{y})^2 + \sum_{i=1}^n (y_{i}-\hat{y}_{i})^2 \end{align} I have a set of exponentially distributed data (x,y) as below,               x   y            332 7.283650     342 7.231924     356 6.949199     369 7.360927     369 7.154024     315 7.379831     334 7.278457     339 7.217902     321 7.238676     300 7.282819     330 7.255710     329 7.126188     374 7.042991     353 7.292686     335 7.405174     360 7.196402     351 7.130031     357 7.218629     348 6.991577     327 7.326131     347 7.453576     335 7.292167     391 7.199726     351 7.310863     307 7.206269     349 7.125773     340 7.129408     341 7.093262     358 7.375157     306 7.516519     301 7.304112     342 7.133147     350 7.345347     333 7.546433     318 7.192559     321 7.142807     347 7.167319     327 7.197025     352 7.122761     343 7.194221     314 10.388952     332 7.425843     344 8.922770     367 7.209697     361 7.040914     309 7.236910     322 7.050262     321 7.974467     291 12.036381     264 10.103215     322 7.208139     313 8.114479     374 7.072074     343 7.376195     333 7.081941     356 7.061375     351 7.176979     349 7.145715     348 6.948888     299 13.043680     337 7.679382     331 7.025437     387 7.182068     349 7.174590     293 7.387724     335 7.240338     330 7.540512     321 9.081374     350 8.787224     270 12.606609     283 8.027543     336 7.113309     346 7.149870     317 10.554307     308 8.318993     317 7.115282     347 7.038421     309 7.059298     344 7.025853     339 7.194532     262 12.690118     291 9.464123     317 12.085718     335 7.581643     348 7.195571     231 8.631216     274 7.765591     298 12.136716     221 8.449553     190 11.299447     285 7.610830     250 9.075558     372 6.836608     312 9.058835     296 6.974543     332 7.381492     298 7.987658     312 7.734847     316 6.975374     227 13.573608     305 9.040347     280 8.783381     362 7.320523     311 6.930607     222 9.994258     299 10.843992     256 9.309050     294 9.304895     303 11.077899     181 14.128775     206 11.550493     293 8.223747     315 8.525272     289 10.906415     279 8.296662     259 9.490090     305 8.418497     267 9.286926     326 8.573466     320 7.391983     20  46.517206     18  46.796088     50  33.885254     10  53.022687     50  33.730700     46  33.914648     30  39.763380     50  34.076680     40  40.189545     66  27.821829     39  40.119954     39  39.077237     57  27.396287     29  46.103090     30  39.334099     29  40.155373     30  40.455963     47  33.197760     60  39.389356     29  40.129821     40  40.059608     38  39.861015     30  46.634160     19  45.749527     47  33.575315     68  33.705149     69  33.255718     69  33.722495     10  52.855150     18  46.515025     28  40.819288     28  39.971840     20  47.167931     28  41.367704     50  34.006985     49  33.528368     50  33.917452     39  39.867870     50  33.280126     40  39.621913     20  46.517206     18  46.796088     50  33.885254     10  53.022687     50  33.730700     46  33.914648     30  39.763380     50  34.076680     40  40.189545     66  27.821829     39  40.119954     39  39.077237     57  27.396287     29  46.103090     30  39.334099     29  40.155373     30  40.455963     47  33.197760     60  39.389356     29  40.129821     40  40.059608     38  39.861015     30  46.634160     19  45.749527     47  33.575315     68  33.705149     69  33.255718     69  33.722495     10  52.855150     18  46.515025     28  40.819288     28  39.971840     20  47.167931     28  41.367704     50  34.006985     49  33.528368     50  33.917452     39  39.867870     50  33.280126     40  39.621913     20  46.517206     18  46.539745     48  33.885254     10  53.022687     50  33.730700     46  33.592453     30  39.763380     50  34.076680     40  40.189545     64  27.515942     36  40.119954     38  39.089597     56  27.260741     28  46.103090     30  39.334099     28  40.155373     30  40.455963     46  33.043622     60  39.389356     28  40.129821     40  40.059608     36  39.724430     30  46.645377     18  45.749527     46  33.307547     66  33.456388     68  33.031573     68  33.714705     10  52.855150     18  46.253489     28  40.595767     28  39.724741     20  47.167931     28  41.117905     50  34.021215     48  33.528368     50  33.917452     38  39.867870     50  33.280126     40  39.621913      I am trying a) to find the MSS, ESS and TSS values.   b) validate if the values I got are correct, satisfying TSS = MSS + ESS ## Approach 1 [Best Fit approach [3], Without converting the exponential data to linear data] \begin{align} y_{i} = y, x_{i} = x \end{align} \begin{align} \overline{y}=23.62464471 \end{align} \begin{align} \hat{y}_{i}=Ae^{Bx_{i}}= 47.3826e^{-0.0055x_{i}} \end{align} where A=exp(a) and B=b. a and b obtained from eq.(3) and eq.(4) from [3]. The equations are reproduced below, \begin{align} a &= (\sum_{i=1}^n ln\ y_{i} \sum_{i=1}^n x_{i}^2 - \sum_{i=1}^n x_{i} \sum_{i=1}^n x_{i}\ ln\ y_{i})\ /\ (n\ \sum_{i=1}^n x_{i}^2 - (\sum_{i=1}^n x_{i})^2) \end{align} \begin{align} b &= (n\ \sum_{i=1}^n x_{i}\ ln\ y_{i} - \sum_{i=1}^n x_{i} \sum_{i=1}^n ln\ y_{i})\ /\ (n\ \sum_{i=1}^n x_{i}^2 - (\sum_{i=1}^n x_{i})^2) \end{align} Calculating the values, \begin{align} MSS = \sum_{i=1}^n (\hat{y}_{i}-\overline{y})^2 = 55961.58373 \end{align} \begin{align} ESS = \sum_{i=1}^n (y_{i}-\hat{y}_{i})^2 = 1739.728157 \end{align} \begin{align} TSS = \sum_{i=1}^n (y_{i}-\overline{y})^2 &= 62029.58951 \end{align} Unfortunately MSS + ESS = 57701.31189 != TSS ## Approach 2 [Best Fit approach [3], Converting the exponential data to linear data] \begin{align} y_{i} = ln y, x_{i} = x \end{align} \begin{align} \overline{y}=2.868328 \end{align} \begin{align} \hat{y}_{i}=a+{bx_{i}}= 3.858255231{-0.0055x_{i}} \end{align} where a and b obtained from eq.(3) and eq.(4) from [3]. The equations are reproduced below, \begin{align} a &= (\sum_{i=1}^n ln\ y_{i} \sum_{i=1}^n x_{i}^2 - \sum_{i=1}^n x_{i} \sum_{i=1}^n x_{i}\ ln\ y_{i})\ /\ (n\ \sum_{i=1}^n x_{i}^2 - (\sum_{i=1}^n x_{i})^2) \end{align} \begin{align} b &= (n\ \sum_{i=1}^n x_{i}\ ln\ y_{i} - \sum_{i=1}^n x_{i} \sum_{i=1}^n ln\ y_{i})\ /\ (n\ \sum_{i=1}^n x_{i}^2 - (\sum_{i=1}^n x_{i})^2) \end{align} Calculating the values, \begin{align} MSS = \sum_{i=1}^n (\hat{y}_{i}-\overline{y})^2 = 150.6993684 \end{align} \begin{align} ESS = \sum_{i=1}^n (y_{i}-\hat{y}_{i})^2 = 3.992417623 \end{align} \begin{align} TSS = \sum_{i=1}^n (y_{i}-\overline{y})^2 &= 154.691786 \end{align} Here we able to validate, MSS + ESS = 154.691786 = TSS ## Approach 3 [Least Square Fit approach [3], Without converting the exponential data to linear data] \begin{align} y_{i} = y, x_{i} = x \end{align} \begin{align} \overline{y}=23.62464471 \end{align} \begin{align} \hat{y}_{i}=Ae^{Bx_{i}}= 48.6062e^{-0.0056x_{i}} \end{align} where A=exp(a) and B=b. a and b obtained from eq.(9) and eq.(10) from [3]. The equations are reproduced below, \begin{align} a &= (\sum_{i=1}^n (x_{i}^2 y_{i}) \sum_{i=1}^n (y_{i}\ ln\ y_{i}) - \sum_{i=1}^n (x_{i} y_{i}) \sum_{i=1}^n (x_{i}y_{i}\ ln\ y_{i}))\ /\ ((\sum_{i=1}^n y_{i}\sum_{i=1}^n (x_{i}^2 y_{i}))-(\sum_{i=1}^n (x_{i} y_{i})) \end{align} \begin{align} b &= (\sum_{i=1}^n y_{i} \sum_{i=1}^n (x_{i} y_{i}\ ln\ y_{i}) - \sum_{i=1}^n (x_{i} y_{i}) \sum_{i=1}^n (y_{i}\ ln\ y_{i}))\ /\ ((\sum_{i=1}^n y_{i}\sum_{i=1}^n (x_{i}^2 y_{i}))-(\sum_{i=1}^n (x_{i} y_{i})) \end{align} Calculating the values, \begin{align} MSS = \sum_{i=1}^n (\hat{y}_{i}-\overline{y})^2 = 59282.06044 \end{align} \begin{align} ESS = \sum_{i=1}^n (y_{i}-\hat{y}_{i})^2 = 1622.208368 \end{align} \begin{align} TSS = \sum_{i=1}^n (y_{i}-\overline{y})^2 &= 62029.58951 \end{align} Unfortunately MSS + ESS = 60904.26881 != TSS ## Approach 4 [Least Square Fit approach[3], Converting the exponential data to linear data] \begin{align} y_{i} = ln y, x_{i} = x \end{align} \begin{align} \overline{y}=2.868328 \end{align} \begin{align} \hat{y}_{i}=a+{bx_{i}}= 3.858255231{-0.0055x_{i}} \end{align} where a and b obtained from eq.(9) and eq.(10) from [3]. The equations are reproduced below, \begin{align} a &= (\sum_{i=1}^n (x_{i}^2 y_{i}) \sum_{i=1}^n (y_{i}\ ln\ y_{i}) - \sum_{i=1}^n (x_{i} y_{i}) \sum_{i=1}^n (x_{i}y_{i}\ ln\ y_{i}))\ /\ ((\sum_{i=1}^n y_{i}\sum_{i=1}^n (x_{i}^2 y_{i}))-(\sum_{i=1}^n (x_{i} y_{i})) \end{align} \begin{align} b &= (\sum_{i=1}^n y_{i} \sum_{i=1}^n (x_{i} y_{i}\ ln\ y_{i}) - \sum_{i=1}^n (x_{i} y_{i}) \sum_{i=1}^n (y_{i}\ ln\ y_{i}))\ /\ ((\sum_{i=1}^n y_{i}\sum_{i=1}^n (x_{i}^2 y_{i}))-(\sum_{i=1}^n (x_{i} y_{i})) \end{align} Calculating the values, \begin{align} MSS = \sum_{i=1}^n (\hat{y}_{i}-\overline{y})^2 = 155.7545382 \end{align} \begin{align} ESS = \sum_{i=1}^n (y_{i}-\hat{y}_{i})^2 = 4.053637091 \end{align} \begin{align} TSS = \sum_{i=1}^n (y_{i}-\overline{y})^2 &= 154.691786 \end{align} Unfortunately, MSS + ESS = 159.8081753 != TSS Questions:   1. Is the above equation is limited to linear data only?   2. How to calculate TSS and ESS for exponentially data without converting it to linear first? The TSS equation seems to be generic that could fit any type of data.  Hope you can guide me on the right way to calculate TSS and ESS for exponential data.               Reference:     [1] http://en.wikipedia.org/wiki/Explained_sum_of_squares     [2] http://www.originlab.com/forum/topic.asp?TOPIC_ID=4823      [3] http://mathworld.wolfram.com/LeastSquaresFittingExponential.html


title: What is the generalized linear model for normal distribution with both mean and variance unknown?

text: When a univariate Gaussian distribution has known variance and unknown mean, its canonical link function is identity. So its generalized linear model under the canonical link function is simply $E(Y|X) = X^T \beta$. When a univariate Gaussian distribution has both mean and variance unknown,   * its canonical link function is complicated (look at "Log-partition $A(\boldsymbol\eta)$" column and "normal distribution" row in the Table) and there are two natural parameters $(\eta_1, \eta_2)$ involved. How is its generalized linear model look like then? I didn't find it on Google, and is it used often?   * Is $Y= X^T \beta + \epsilon, \epsilon \sim N(0, \sigma^2)$ with unknown parameter $(\beta, \sigma^2)$ still a generalized linear model?  Thanks!


title: How to test for mean and variance both affected by a variable?

text: Suppose you observe $\left(x_i,y_i\right)$, where $y_i \sim \mathcal{N}\left(\mu_i,\sigma_i^2\right)$. There are a number of ways that the $\mu_i, \sigma_i$ might vary with $x_i$:   1. If one suspected $\mu_i$ varied with $x_i$, but $\sigma_i$ is independent of $x_i$, one could perform plain old linear regression.   2. If one suspected $\sigma_i$ varied with $x_i$, a heteroskedasticity test ( _e.g._ Breusch-Pagan) could be performed. I am interested in a third case where ${\mu_i}/{\sigma_i}$ vary linearly with $x_i$, but there is no correlation assumed between $\mu_i$ and $x_i$, and no correlation assumed between $\sigma_i$ and $x_i$. Is there a standard test for this problem?


title: How do you select the best time series model?

text: If you are given the parameter estimates and their p value and ACF PACF only how to can you chose the most suitable model if there are more than one model where parameters are significant?


title: Assessing the significance of differences in distributions

text: I have two groups of data. Each with a different distribution of multiple variables. I'm trying to determine if these two groups' distributions are different in a statistically significant way. I have the data in both raw form and binned up in easier to deal with discrete categories with frequency counts in each. What tests/procedures/methods should I use to determine whether or not these two groups are significantly different and how do I do that in SAS or R (or Orange)?


title: Why is dependence a problem?

text: I'm interested in why dependent observations are a problem in statistics. Let's say you want to know if there is a difference in mean exam scores between two schools. You collect 50 observations in each school. These 50 observations are derived from 5 different classrooms in each school and there is dependence within classrooms. In this instance, how would the results of the t-test be affected and how might they lead to inaccurate conclusions?


title: Using graph to depict individual movement pre and post intervention

text: I had an intervention with control group and experimental group. Both groups had pre- and post-tests. So I was using 4 graphs to explain in a presentation. (pre-test control, pre-test experiment & post-test control and post-test experiment) Each graph showed responses to a likert scale question with 5 options. One of the listeners asked me if I could show the shift of 'how many' students had moved from one option in pre-test to some other option in post-test. I can use pivot tables in Excel however, I am not very good in stats. Can someone kindly guide me through this?


title: Whether to use Spearman's rho or multiple regression to examine relationship between two Likert scales?

text: _I am doing my thesis and have absolutely no previous experience in statistics._ I have constructed several Likert scales by forming composites scores each based on 4-6 items which tests the level of agreement of my respondents. Specifically, I have created two scales 'purchase behaviour' and 'website appeal'. I want to see whether 'website appeal' is correlated with 'purchase behavior'. Now that I am doing my analysis, I am confused as to whether I should use Spearman rho or multiple regression. The Spearman correlation and multiple regression have different p-values, so much so that one states that I need to reject my hypothesis and another accept. **So in this case should I use Spearman's rho or multiple regression?** Is there a theoretical rule that I must use, say, multiple regression because I am testing 4-6 items on the likert, although I have grouped them together and intend to 'view' them as two single variables. * * * Thanks, Jeromy. The article by Gelman and Stern (2006) is really interesting! Being the non-statistician me and trying to get on with my MBA thesis, I would be very tempted to find an analysis which gives me a simple method to analyse my data and ultimately, test my hypotheses. I know this shouldn't be the way, but stats aren't exactly fun nor interesting. I was talking to my supervisor and he suggested using regression when I had planned to use Spearman (cos Likert scale items are ordinal and if I want to test ordinal vs ordinal, I use Spearman - according to my research methods text) Yes, I am planning to test only two variables (predictors?) at a time, so, Spearman can technically be used. But these two variables (both dependent and independent) are computed as a new variable from the different items I have on the Likert scale (does this make sense?). I am just concerned that my analysis would be deemed incorrect if I used the 'wrong' statistical analysis - or doesn't this matter?


title: What does R do when it plots the residuals of an AR fit?

text: This is a question that's been bugging me for some time. The problem is this: I'm modelling the residuals of a model $f(t,\vec{\theta})$ with (what I think is) an AR process plus a white noise process via MCMC using a multivariate gaussian likelihood, where I model the covariance matrix elements equal to the elements of the autocovariance of an $AR(1)$ process plus i.i.d. white noise $\varepsilon(t)$, i.e., I model $$r(t)=AR(1)+\varepsilon(t),$$ where $r(t)$ are the residuals of my model, $r(t)=d(t)-f(t,\vec{\theta})$, where $d(t)$ is the data. The thing is that I wanted a visual check to see if my fit was ok (because I'm pretty sure via some previous analysis that my residuals can be modelled efficiently by an $AR(1)$ plus white noise: I just want to show this result to other people "quickly"). What I did then was to analize the residuals via the `arima` and `arma` functions in the `tseries` library in R, and I saw that they actually fitted the residuals with the same coefficient as the MCMC result (within the error bounds and not taking in consideration the additive white noise, which has small variance anyways). However, I also saw that the `arma` function has some nice plots: it plots the "residuals" of the AR fit. This is my question: **how do you plot the residuals of an AR fit if an AR model is stochastic by nature?** I think my question goes down to: **how do you generate a realization of your model that matches your particular realization of the AR process (i.e. your data) in order to substract this realization from your data?** The only way I could think of is to simulate various realizations of the $AR(1)$ model until one fits your data, but this seems like stacking the deck.


title: Looking for example of non-elliptical, multivariate, absolutely continuous distributions with non independent components

text: I would like to test several (multivariate --between 5 and 20 dimensions) depth measures empirically. I'm looking for a couple of test-benches (read distributions). These should be:   1. non-elliptical,    2. absolutely continuous,   3. uni-modal    4. have convex contours (if you prefer, level sets)   5. have non independent components.    6. Easy to implement (i.e. implemented in R or source code available on-line).    7. defined in arbitrary dimensions  Any suggestions?


title: Reference for hierarchical Bayesian modelling

text: I am currently reading "Bayesian Data Analysis" by Gelman et al. and my main goal was to learn about Hierarchical modelling on chapter 5. I read until chapter 4 and the book is written terribly for a taste of a math student as it is pretty sketchy and engineering oriented. I decided to not to continue anymore with this book and I would be very grateful if somebody could introduce a reference with a more rigorous approach to the topic.


title: Time Varying System Matrices in Kalman Filter

text: Kalman filter can accommodate time varying system matrices. Equations to run the filter are the same and it preserves its optimality under linear gaussian model. My question is the following: Can the evolution of time varying system matrices be stochastic? In some references I seem to read between the lines that they should evolve deterministically. Does it mean that the entire filter breaks or do we simply lose optimality by making them stochastic? For reference, please peek at section 3.2 of the following paper: http://www.ims.cuhk.edu.hk/~cis/2012.1/CIS%2012-1-05.pdf A similar comment is in Harvey's book on Kalman Filter. Thanks


title: Log Likelihood for GLM

text: In the following code I perform a logistic regression on grouped data using glm and "by hand" using mle2. Why does the logLik function in R give me a log likelihood logLik(fit.glm)=-2.336 that is different than the one logLik(fit.ml)=-5.514 I get by hand?               library(bbmle)          #successes in first column, failures in second     Y <- matrix(c(1,2,4,3,2,0),3,2)          #predictor     X <- c(0,1,2)          #use glm     fit.glm <- glm(Y ~ X,family=binomial (link=logit))     summary(fit.glm)          #use mle2     invlogit <- function(x) { exp(x) / (1+exp(x))}     nloglike <- function(a,b) {       L <- 0       for (i in 1:n){          L <- L + sum(y[i,1]*log(invlogit(a+b*x[i])) +                     y[i,2]*log(1-invlogit(a+b*x[i])))       }      return(-L)      }            fit.ml <- mle2(nloglike,                start=list(                  a=-1.5,                  b=2),                data=list(                  x=X,                  y=Y,                  n=length(X)),                method="Nelder-Mead",                skip.hessian=FALSE)     summary(fit.ml)          #log likelihoods     logLik(fit.glm)     logLik(fit.ml)               y <- Y     x <- X     n <- length(x)     nloglike(coef(fit.glm)[1],coef(fit.glm)[2])     nloglike(coef(fit.ml)[1],coef(fit.ml)[2])


title: Generalized Linear Mixed Models: Diagnostics

text: I have a random intercept logistic regression (due to repeated measurements) and I would like to do some diagnostics, specifically concerning outliers and influential observations. I looked at residuals to see if there are observations that stand out. But I would also like to look at something like Cook's distance or DFFITS. Hosmer and Lemeshow (2000) say that due to the lack of model diagnostic tools for correlated data, one should just fit a regular logistic regression model ignoring the correlation and use the diagnostics tools available for regular logistic regression. They argue that this would be better than doing no diagnostics at all. The book is from 2000 and I wonder if there are methods available now for model diagnostics with mixed effects logistic regression? What would be a good approach to check for outliers? **Edit (Nov 5, 2013):** Due to the lack of responses, I am wondering if doing diagnostics with mixed models is not done in general or rather not an important step when modeling data. So let me rephrase my question: What do you do once you found a "good" regression model?


title: Classification tree with rpart package

text: I am using R with package `tree` to build a classification tree to summarize and represent my data. I have only 27 records and 10 variables. I built the tree but the but the `xerror` grows:                  CP  nsplit  rel error  xerror     xstd     1 0.2       0          1     1.0  0.40369     2 0.0       5          0     1.4  0.45542      Is due to overfitting problem? unbalanced data? Can I still use the tree to describe the data?


title: Longitudinal design - repeated measures, non-parametric, non-balanced design

text: I am currently trying to find the best way to analyse my data. Participants were tested longitudinally on up to five occasions. On each day participants were assessed for the presence or absence of a neurological marker (dichotomous), and were tested on a continuous measure (scoring 0-10). I simply want to compare how the positive or negative diagnosis on any given occasion predicted score on the continuous measure. Unfortunately as the study took place in a critical-care hospital ward participants were all tested a different number of times, and were positive or negative for varying proportions of these observations. The distributions of scores for each group on the continuous variable were also non-normal. So, I was wondering if there was a non-parametric repeated measures test that could handle unbalanced designs? I have been told Multi-level modelling might be an option, but I don't have much experience with it. Thanks in advance for any help!


title: Binary classification problem

text: I have a binary classification problem. My inputs consist of a time-series of values plus some binary values. For real-valued inputs I would usually use a neural network, while for binary-values inputs I would use either a neural network with a step function or some other method fit for binary inputs (like a maximum entropy classifier) In the case where input types are mixed what architecture should be used?


title: Post-hoc analysis for Logrank test

text: I need to perform survival analysis for 4 different groups in a study. When I use the Logrank test, the null hypothesis is rejected. My question is, how can I tell now exactly which pair of groups causes the null hypothesis to be rejected? Is there a special post-hoc test, or should I repeat the logrank test manually for different pairs of groups? (or another solution, if possible) Appreciate any kind of help!


title: How would one calculate intrarater reliability in Stata?

text: My question follows from a previous excellent response in another post: Intra- and inter-rater reliability on the same data Wow, what a truly great response! I am running through a similar problem and I wondered if someone would be able to expand this example to a dataset similar to the one you created **but** with 2 measurements per patient, per rater, to allow for calculation of intra-rater reliability? And with scoring based on a 1-6 scale instead of a binary response? The dataset would look like this:               +-------------------------------------------------------------+     | a1   a2   b1   b2   c1   c2   d1   d2   Time  Patient |     |-------------------------------------------------------|     |  3    2    1    2    1    4    2    1     1       1   |     |  4    4    3    5    1    5    2    1     2       1   |     |  1    2    2    5    3    5    1    3     1       2   |     |  2    1    3    5    2    6    1    1     2       2   |     |                                                       |     +-------------------------------------------------------------+


title: Difference in proportions but one of them is 0

text: What is the best way to test for difference of proportions when of them is 0% (or 100%) ? edit: I have two raters, who give a score of either 1 or 0. One of them gave a 0 in every single case. Now I want to know if the proportions of 1's is significantly different between the two.


title: Sampling the outcome of a function in R

text: I have a function which simulates the outcome of a particular card game, and returns a value denoting the outcome of the game. I want a concise syntax for sampling multiple plays of the game, and getting a vector containing the outcomes. This is what I want to be able to write:               sample(cardgame(params), iterations, replace=?)      which would ideally return a vector like (0, 1, 1, 1, 0, 1, 0, ...). The above syntax actually evaluates cardgame(params) once, and treats it as numerical parameter to sample. Is there a way to avoid creating the outcomes vector by hand with a for loop and multiple calls of the cardgame function? I ask, because being able to call sample() on user-created distributions would be very powerful. Thanks!


title: So the interpretation of 2-way interaction cannot be interpreted if it is significant in 3-way interaction?

text: I use LMM to analyse my data and my variables are as follows: DV: continuous IV: 1) color(red, blue, green), height(low, tall, medium), and sex (male, female) Then for color, "red" is baseline, and for height, "low" is baseline. The output is like this:               colorblue (sig)     colorgreen     heighttall     heightmedium     male     colorblue:heighttall (sig)     colorgreen:heighttall     colorblue:heightmedium     colorgreen:heightmedium     colorblue:male     colorgreen:male     heighttall:male     heightmedium:male     colorblue:heighttall:male (sig)     colorgreen:heighttall:male     colorblue:heightmedium:male (sig)     colorgreen:heightmedium:male (sig)      For the interpretation, is it appropriate to say: The colorblue has a main effect when red=0, and the differences between colorblue and colorred is lower (or higher depending on coefficient) in heighttall compared to heightlow. In 3-way interaction, the effect of colorblue and heightall is greater (or lower depending on coefficient) in male compared to female. However, as the interaction of colorblue:heightall:male is significant, this implies that the effect of color alone or color and height is insufficient. In conclusion, the effect of color correlates to height and sex. Is my interpretation correct?


title: Effect size for interaction effect in pre-post treatment-control design

text: If you choose to analyse a pre-post treatment-control design with a continuous dependent variable using a mixed ANOVA, there are various ways of quantifying the effect of being in the treatment group. The interaction effect is one main option. In general, I particularly like Cohen's d type measures (i.e., ${\frac{\mu_1 - \mu_2}{\sigma}}$). I don't like variance explained measures because results vary based on irrelevant factors such as relative sample sizes of groups. Thus, I was thinking I could quantify the effect as follows   * $\Delta\mu_c = \mu_{c2} - \mu_{c1}$    * $\Delta\mu_t = \mu_{t2} - \mu_{t1}$    * Thus, the effect size could be defined as $\frac{\Delta\mu_t - \Delta\mu_c}{\sigma}$  where $c$ refers to control, $t$ to treatment, and 1 and 2 to pre and post respectively. $\sigma$ could be the pooled standard deviation at time 1. **Questions:**   * Is it appropriate to label this effect size measure `d`?   * Does this approach seem reasonable?   * What is standard practice for effect size measures for such designs?


title: Pooling samples -- what are the effects

text: There is an experiment in which biological samples are taken for several treatments and a control. For each treatment, we have enough resources to evaluate 5 samples. Some researchers suggest (for technical reasons) that for each treatment, we should collect ten samples and pool them, such that we have five pools, each consisting of two biological samples mixed together. My main interest is in understanding the treatment effect size in the context of biological variability between the samples, which I think will be lost (or complex to analyse) if each pool will be an average of two biological samples (there are also other reasons for not pooling, mostly due to the fact that the pools will be incompatible with other samples and measurements which are not to be pooled). However, the possibility exists that due to technical reasons we will be forced to pool. **My questions:**   1. Can we estimate the underlying biological variance from the variance of the pools? I think that since we the expected variance of the pools will be $\frac{\sigma^2}{2}$ (given that real the biological variance in the population $\sigma^2$), could we multiply the variance estimate from the pools by $2$ to get the estimate of the biological variance, correct?   2. The experimental design is actually a bit more complex and it's analysis will involve modelling group and random effects (a mixed effect model). The main question will be to understand the treatment effects in the context of the population (and not of the pools). Is it possible to incorporate the pooling into the model? To elaborate on the second question, I will give a specific example of a model. Consider a randomized block design as follows (I do hope I get it right...): $Y_{ij}=\beta_0 + \beta_j\cdot T_j + \eta_i\cdot G_i + \epsilon_{ij} $ in which there is a fixed effect of the treatment ($\beta_j$) and a random effect of the experimental batch ($\eta_i$). I would analyse it using the `nlme` or `lmer` libraries from R, applying mixed effect models. But I am confused as to what I would have to do if the samples were pooled.


title: VC dimension of SVM with polynomial kernel in $\mathbb{R^{2}}$

text: What is the VC dimension of SVM with the polynomial kernel $k(x,x')=(1+<x,x'>_{\mathbb{R^{2}}})^{2}$ for binary classification in $\mathbb{R^{2}}$? It would be equal or more than v iff there exists a set of v points such that, any labeling (-1 or +1) of the points given, there exists a correct separating border. It would be strictly less than v iff for all sets of v points, there exists a labeling of the points such that there is no correct separating border. In this case, the separating border is a conic section or a line, so any idea based on these curves rather than SVM is welcome. ~~For instance, let $x=(x_{1},x_{2}) \in \mathbb{R^{2}}$. It is mapped to $\mathbb{R^{6}}$ using a certain $\phi$ and an hyperplane in this new space is a conic section in the former space. So my guess would be the VC dimension of a linear classifier in $\mathbb{R^{6}}$, which is 7.~~ * * * Actually, **the answer is 6** : the equation of an hyperplane in $\mathbb{R^{d}}$ is $$<w,x>_{\mathbb{R^{d}}}+b=0$$ where $w \in \mathbb{R^{d}}$ and $b \in \mathbb{R}$. Here, we have $b=1$ and $w \in \mathbb{R^{5}}$, not $\mathbb{R^{6}}$.


title: Interpreting Two-way repeated measures ANOVA results: Post-hoc tests allowed without significant interaction?

text: I have conducted a study where I employed a two-way repeated measures design to investigate whether subjects respond to a main treatment effect. What I'm mainly interested in is whether subjects subsequently compensate for the treatment effect by returning to a control level, or whether the initial differences persist for the duration of the treatment (time x treatment interaction). My factors: Treatment: Between-subjects factor (3 levels) Exposure Time: Within-subjects repeated measure (4 levels - pre-exposure values not included) Treatment x Time Interaction: ![Response over time to treatment effect \(3 levels\). Pre-test \(time zero\) values were not included in the ANOVA](http://i.stack.imgur.com/w867m.jpg) Error bars = Standard deviation My ANOVA results: Sphericity was taken into account: p = 0.25. Adjusted values told the same story. Both main effects, 'Treatment' and 'Time' were significant. Treatment effect had a large effect size (0.4). (sensu Cohen) Time, isn't really that interesting on its own, had a small effect size (0.08) The interaction was not significant (p = 0.3) and had a tiny effect size (0.04). My problem: The graph appears to show that one of the treatment levels interacts with exposure time (fig. d). (ie the lines are not all parallel from day 4 - 14) The ANOVA results however suggest that that the difference between treatments persists (the lines remain parallel). I have looked at pairwise comparisons within levels of the main effect. My question: I am tempted to included post-hoc comparisons to compare treatment levels at each time interval to show that the difference between groups disappears on day 11. Does this defeat the purpose of running an ANOVA, which clearly tells me that I should just state that there was a main effect (ie the treatment groups remained different in fig d). I therefore shouldn't really fish around for significant differences. The literature that I'm comparing my work to employs ANOVA, but then also puts pairwise comparisons all over the place. I'm a fan of effect sizes, rather than p-values, so feel justified to just highlight the main treatment effect. My main conclusion could be very different depending on whether I go with the ANOVA (ie no interaction -> no recovery) or whether I include pairwise comparisons at each time point (partial recovery). Side note: Mixed effects models in lme gave the same result. I stuck with ANOVA as my audience is more likely to relate to it. Thanks for the feedback.


title: Treating the log-likelihood as a probability distribution and normalizing?

text: I have modeled a distribution, $f$, over a r.v. $x \in \mathbb{R}^3$. At inference a set of measuring points, $X$, of the r.v. variables show up. I want to form a distribution over this sample set so that I can sample from it so that more probable points will show up more often. My idea was to evaluate $f$ for each point in the set and then normalize it to form a discrete distribution over the points. Since $f \approx 0$ I want to use the log likelihood. However, the $\log$ is a nonlinear transform and normalizing it means calculating: $\log \sum_i f(x_i)$, which is bound to give numerical errors due to the fact that $f(x_i)\approx 0$. So the best I can do is a log sum exp calculation. My question is: is there a better way of doing this? And is this a completely wrong approach in terms of modeling probability distributions?


title: Classification with Bayes Network (deal)

text: I am working with multiple binary vectors e.g., A,B,C,D,E,F,G,H. I want to find the classification between them. I have tried the following:               log_data<-read.csv(choose.files(), as.is = T, header = T, blank.lines.skip = TRUE)     data<-log_data[2:ncol(log_data)]     data     TIME  A B C D E F G     1     1 1 1 0 1 0 1 1     2     0 0 1 1 1 1 0 1     3     1 1 1 1 1 0 1 1     4     1 0 1 1 1 1 0 1     .....................     fit <- network(data)     fit.prior <- jointprior(fit)     fit <- getnetwork(learn(fit,rats,fit.prior))     Error in postc0c(node$condposterior[[1]]$mu, node$condposterior[[1]]$tau,  :      NA/NaN/Inf in foreign function call (arg 1)      Getting this error just because all are continuous variable and NULL at mu. How should I proceed in order to classify after creating a network?


title: How to extract the weight and bias in libsvm?

text: I want to extract the weight and the bias for my dataset which has two columns and more than 3000 rows. Also I want to draw a line by this equation: y = w.x + b where w is the weight, x is the data point, and b is the bias. Could anyone help me, please? Regards,


title: Fisher's exact test, very small sample size, RxC table, gave me bad result

text: I apply Fisher's exact test on very very small sample size. Each location in the table would have about 1 or two samples. The fisher's method has no room to enumerate the the tables, and returns me a very bad result. So is there any other improved exact test for very tiny sample size data?


title: Testing goodness of fit

text: I am using MATLAB 2012a to run 3 tests. With my data, I obtain very small values, and for the Anderson–Darling (A–D) tests, the results are not even normalized. Please verify my code. I am also using the A–D script which is found at this link.               sample= N;          h1=histfit(sample,30,'weibull');     xdata1 = get(h1(2), 'XData');     ydata1 = get(h1(2), 'YData');     [~,p_weibull_cs,stats_weibull_cs] = chi2gof(ydata1);     [~,p_weibull_ks,stats_weibull_ks] = kstest(ydata1);     p_weibull_ad = AnDartest(ydata1);


title: Test whether there is a significant difference between two groups

text: I have some trouble with the last test I need to perform for my bachelor thesis. I have two variables and I need to test whether there is a significant difference. The first variable has all the returns of stocks on specific days. The second variable has these returns as well, however some days have been deleted. I want to test whether there is a significant difference between 'before' and 'after' the filtering. Is it OK if I perform a Paired Sample T-test, or should I use another test? Every help is very appreciated!


title: Pearson correlation coefficient using MATLAB when standard deviation of one variable is zero

text: For calculating Pearson Coefficient, i am using the below formula in MATLAB:               C=cov(x,y);     p=C(2)/(std(x)*std(y));      But C(2) is zero in my case if               x = [1 0 1 0 0 0];     y = [1 1 1 1 1 1];      Any idea where I am going wrong


title: Under what conditions will Kullback-Leibler divergence/mutual information be infinity?

text: For two perfectly correlated Gaussian variables, the mutual information between them, and thus the KL divergence between the product of the marginal distributions and the joint distribution, is infinity. I have read that KL divergence is always finite if Q is continuous in respect to P, but I am having a hard time understanding why, and how to know what the maximum KL divergence between the product of the marginals and the joint is for any given marginal distributions. I have a feeling that the maximum possible mutual information between two distributions is related to the KL divergence of their marginal distributions - i.e. a uniformly distributed variable cannot have infinite mutual information with a Gaussian distributed variable, and that is because the KL divergence between two two marginals is not 0. However, if I know the KL divergence of two marginal distributions, is it possible to calculate the maximum possible mutual information between them? I hope this isn't a stupid question, but I know there must be a relationship and I can't find anything about it readily. Thanks!