title: Is this a case for an ordinal logistic regression? Problems interpreting output

text: I'm a beginner in statistics and R, sorry if this question may seem trivial. I've collected data measuring several different parameters in 40 subjects at two time-points (t1 and t2). There are 3 main parameters in which I'm interested, let's call them ParA, ParB, ParC. ParA is a score of disability. It is on an arbitrary scale (so it is an ordinal scale measure, if my understanding is correct) and values range from 0.0 to 10.0. Note that the increments in this scale are by 0.5 unit, so values like, e.g. 1.5 are possible. I have two measures, at t1 and t2, so I can describe at least three variables from ParA: ParA at t1, ParA at t2, and whether a subject progressed or not (0 or 1). Being a ratio scale measure, I think it would not make much sense to compute a difference (eg. ParA at t2 - ParA at t1), but I'm willing to accept suggestions on this matter. ParB and ParC are meausurements of two anatomical structures; they are continuous; Par B is an area measured in mm2, ParC is a volume measured in mm3. I believe they should be considered ratio scale measure. For each one I can describe at least 4 variables: measuments at t1 and t2 (eg. ParB1, ParB2), absolute difference between the two measurements (ParB2-ParB1), percentual difference between the two. What I want to do: I want to do some sort of regression, to see if ParA at t2 is best predicted by B or C (percentual difference, or maybe ParB or ParC at t1, I've yet to decide which makes most sense), but I don't know how to do this. I tried doing some research and concluded that I want to do an ordinal logistic regression, but now I'm unsure how to interpret the results and I'm questioning my choice. Here's some data in case someone wants to reproduce my situation:               require(MASS)          example.df <- data.frame(ParA1=c(1.5,0.0,0.0,1.0,1.5,1.0,1.0,0.0,0.0,0.0,3.0), ParA2 = c(2.5,1.5,1.0,2.0,2.0,1.5,2.0,0.0,0.0,0.0,6.5), Progressed=c(1,1,1,1,1,1,1,0,0,0,1), ParB1=c(222.76,743.07,559.65,642.93,584.36,565.53,590.88,465.31,570.22,543.91,574.80), ParB2=c(214.5,674.71,538.75,560.72,581.9,566.40,499.72, 434.72,528.33,517.61,516.1), ParBAbsolDiff=c(-8.27,-68.36,-20.90,-82.21,-2.46,0.87,-91.16,-30.59,-41.88,-26.31,-58.71), ParBPercentDiff=c(-3.71,-9.20,-3.73,-12.79,-0.42,0.15,-15.43,-6.57,-7.34,-4.84,-10.21), ParC1=c(1585354,1600993,1818728,1595059,1445126,1599984,1454398,1540987,1567783,1559505,1523271), ParC2=c(1578834,1512068,1800791,1514774,1472185,1548337,1440284,1505046,1586734,1622379,1496734), ParCAbsolutDiff=c(-6520.26,-88925.62,-17937.04,-80285.40,27059.77,-51646.81,-14114.52,-35940.91,18951.04,62873.71,-26536.51), ParCPercentDiff=c(-0.41,-5.55,-0.99,-5.03,1.87,-3.23,-0.97,-2.33,1.21,4.03,-1.74))          > myregression <- polr(ParA2 ~ ParBPercentDiff + ParCPercentDiff,data=example.df, Hess=TRUE)     Error in polr(ParA2 ~ ParBPercentDiff + ParCPercentDiff, data = example.df,  :      response must be a factor          > example.df$ParA2 <- factor(example.df$ParA2)     > myregression <- polr(ParA2 ~ ParBPercentDiff + ParCPercentDiff, data=example.df, Hess=TRUE)          > summary(myregression)     Call:     polr(formula = ParA2 ~ ParBPercentDiff + ParCPercentDiff, data = example.df,      Hess = TRUE)          Coefficients:                    Value Std. Error t value     ParBPercentDiff -0.04825     0.1114 -0.4330     ParCPercentDiff -0.13650     0.2079 -0.6566          Intercepts:         Value   Std. Error t value     0|1     -0.4982  0.9546    -0.5219     1|1.5   -0.0267  0.9367    -0.0285     1.5|2    0.7736  0.9874     0.7835     2|2.5    2.1062  1.1628     1.8113     2.5|6.5  2.8957  1.3531     2.1400          Residual Deviance: 35.89846      AIC: 49.89846      > ci <- confint(myregression)     Waiting for profiling to be done...     > exp(cbind(OR= coef(myregression), ci))                        OR     2.5 %   97.5 %     ParBPercentDiff 0.9528960 0.7596362 1.200121     ParCPercentDiff 0.8724038 0.5670611 1.321134      I searched the internet, but I don't understand what I found. There are several questions in stackoverflow, but the answers were too advanced. I'm a medical student, so I have some basics concept in statistics, but this is advanced stuff for my current level of understanding. I'm trying to study statistics in my free time, but at the moment I'm pressed, and I need to understand what this mean; I think it is best if it is explained in simple terms. Getting back to my question... I'm confused by the results I obtained, probably because I'm confused about some underlying concept. What I conclude from looking at summary(myregression) is that ParBPercentDiff is a worse predictor than ParCPercentDiff (a certain increase in ParBPercentDiff gives a decrease of -0.04 in the expected value of ParA2 on the log odds scale - so a decrease of ParBPercentDiff should give an increase of +0.04; these values are higher for ParCPercentDiff); however the OR seems to tell a different story; Like, the odds of ParA2 increasing are greater with an increase in ParBPercentDiff...? My interpretation is surely flawed. I think I could interpret the results more easily if I could plot this regression, but I didn't get what I expected.               myprof <- profile(myregression)     plot(myprof)      So:   * Is ordinal logistic regression really what I want to do in this case?   * If not, what analysis do you suggest?   * If yes, how can I interpret the result I got (please, in simple terms)?


title: Conjugate Prior for Probit likelihood function

text: I am trying to do a Bayesian analysis in which my likelihood function is a probit function on two parameters. From various sources, I found out that Normal distribution is a conjugate prior to probit likelihood, but I am not being able to compute the posterior from the functions. Here are the details of the problem Likelihood: $\Phi(\theta-\beta)$ Prior for $\beta$: $\phi(\beta-\mu_1)$ Prior for $\theta$: $\phi(\theta-\mu_2)$ Joint Posterior: $C\Phi(\theta-\beta)\phi(\beta-\mu_1)\phi(\theta-\mu_2)$ Now I want to get the marginal distributions from this joint posterior


title: AIC BIC Mallows Cp Cross Validation Model Selection

text: If you have several linear models, say model1, model2 and model3, how would you cross-validate it to pick the best model? (In R) I'm wondering this because my AIC and BIC for each model are not helping me determine a good model. Here are the results:               Model - Size (including response) - Mallows Cp - AIC - BIC      Intercept only - 1 - 2860.15 - 2101.61 - 2205.77     1 - 5 - 245.51 - 1482.14 - 1502.97      2 - 6 - 231.10 - 1472.88 - 1497.87       3 - 7 - 179.76 - 1436.29 - 1465.45       4 - 8 - 161.05 - 1422.43 - 1455.75       5 - 9 - 85.77 - 1360.06 - 1397.55       6 - 10 - 79.67 - 1354.79 - 1396.44      7 - 17 - 27.00 - 1304.23 - 1375.04     All Variables - 25 - 37.92 - 1314.94 - 1419.07      Note - assume the models are nested.


title: How do I tell when and why to use specific statistical measures?

text: I've taken a few probability classes and now understand how to calculate some statistical measures like mean and confidence intervals. What I don't know is the what, when, and why of using these measures for specific situations. I'm hoping to put together a good collection of each of these measures, what they're used for, and what situations these are good to use. Specifically I'm looking for these (but not limited to):   * Mean (average)   * Standard Deviation   * Variance   * Confidence Intervals   * Median


title: How to calculate SE for a binary measure, given sample size n and known population mean?

text: Lets say there is population of measurements X, and 50% of those X = 1, and the the other 50% = 0. Therefore, population mean = 0.5 given a random sample of size n, how do you determine the SE of that sample? Or, in layman's terms, if you flip a coin n times, how far can you expect to deviate from 50% heads and 50% tails?


title: How to fit Gaussian

text: I have a series of (x1,y1) points how to get back a series of (x2,y2) points that are on a Gaussian curve that has been fit to the data.


title: Data transposition from 'clustered rows' into columns

text: I am facing a difficult challenge, given my very low skills at text mining… Basically, I have a list of approx. 200 individuals described in a plain text file following a simple structure: > `N: (name)`   >  `Y: (year of birth)`   >  `S: (sibling)` > > `N: (name)`   >  `Y: (year of birth)`   >  `S: (sibling)` > > `[etc.]` Here's the catch:   * each individual can have missing data, in which case the line is _omitted_ ;   * multiple siblings, in which case his/her entry has more than one `S:` line. The only real constants in the text file are:   * each individual is separated by a blank file from each other;   * all lines feature the appropriate `N:`, `Y:` or `S:` prefix descriptor. Is there a method, using pretty much anything from Excel (!) to Stata to R or even Google Refine or Wrangler, to turn that organised chaos into a standard dataset, by which I mean, one column per descriptor, with `S1, S2... Sn` columns for siblings? If that is the wrong forum, please redirect me. I just figured that the stats crew would be the most acquainted with text mining.


title: Testing for stability in a time-series

text: > Is there a standard (or best) method for testing when a given time-series > has stabilized? * * * ### Some motivation I have a stochastic dynamic system that outputs a value $x_t$ at each time step $t \in \mathbb{N}$. This system has some transient behavior until time step $t^*$ and then stabilizes around some mean value $x^*$ with some error. None of $t^*$, $x^*$, or the error are known to me. I am willing to make some assumptions (like Gaussian error around $x^*$ for instance) but the less a priori assumptions I need, the better. The only thing I know for sure, is that there is only one stable point that the system converges towards, and the fluctuations around the stable point are much smaller that the fluctuations during the transient period. The process is also monotonic-ish, I can assume that $x_0$ starts near $0$ and climbs towards $x^*$ (maybe overshooting by a bit before stabilizing around $x^*$). The $x_t$ data will be coming from a simulation, and I need the stability test as a stopping condition for my simulation (since I am only interested in the transient period). ### Precise question Given only access to the time value $x_0 ... x_T$ for some finite $T$, is there a method to say with reasonable accuracy that the stochastic dynamic system has stabilized around some point $x^*$? Bonus points if the test also returns $x^*$, $t^*$, and the error around $x^*$. However, this is not essential since there are simple ways to figure this out after the simulation has finished. * * * ### Naive approach The naive approach that first pops into my mind (which I have seen used as win conditions for some neural networks, for instance) is to pick to parameters $T$ and $E$, then if for the last $T$ timesteps there are not two points $x$ and $x'$ such that $x' \- x > E$ then we conclude we have stabilized. This approach is easy, but not very rigorous. It also forces me to guess at what good values of $T$ and $E$ should be. It seems like there should be a better approach that looks back at some number of steps in the past (or maybe somehow discounts old data), calculates the standard error from this data, and then tests if for some other numbers of steps (or another discounting scheme) the time-series has not been outside this error range. I included such a slightly less naive but still simple strategy as an answer. * * * Any help, or references to standard techniques are appreciated. ### Notes I also cross posted this question as-is to MetaOptimize and in a more simulation-flavored description to Computational Science.


title: Optimal weighting matrix instrumental variables estimator

text: The formula for the optimal weighting matrix when you perform regression with more instrumental variables than endogenous predictors is the following: $W_{opt} = (\frac{1}{N}Z'Z)^{-1} $ This tells us that we only have to look at the variance covariance matrix of the instruments, But doesn't it make more sense to give more weight to the strongest instruments (or in other words, those who correlate better with the endogenous predictors)? Thanks in advance!


title: Need subspace partition algorithm, not necessarily a full classifier

text: ![enter image description here](http://i.stack.imgur.com/qmNwR.png) The image above represents a hypothetical data set of interest. For some set of points in N-dimensional space (each attribute of the data set corresponds to one dimension), I want to identify partitions with relatively higher density of hits (red points) vs. misses (black points). My goal is not to classify all points. For example, in the data set represented above, there is a partition defined by y > 0.7 and y < 0.75, z > 0.7 that contains a relatively high density of hits. I am looking for an algorithm that can identify such a partition in N-dimensional space (N might range from 5-20). Does anyone know of a good way to do that? If there is a better place to ask this question, please let me know.


title: What is the best way to determine if pageviews are trending upward or downward?

text: Given the following dataset for a single article on my site:               Article 1     2/1/2010 100     2/2/2010 80     2/3/2010 60          Article 2     2/1/2010 20000     2/2/2010 25000     2/3/2010 23000      where column 1 is the date and column 2 is the number of pageviews for an article. What is a basic acceleration calculation that can be done to determine if this article is trending upwards or downwards for days I have pageview data for it? For example looking at the numbers I can see Article1 is trending downwards. How can that be reflected in an algorithm most easily? thanks!


title: Convex loss optimization and supervised learning

text: I once heard the following statement from the web > Convex loss functions and nonlinear classification are two important > concepts in supervised learning. I do not know why this statement can be made? Or what are the roles that convex loss functions and nonlinear classifications play in developing supervised learning algorithms?


title: Burn-in period for random walk

text: We are trying to make simulation experiment involving a common stochastic trend, that is described by the random walk (or $I(1)$ process) $Y_t = Y_{t-1} + \varepsilon_t$, where innovations $\varepsilon_t$ ~ $N(0,1)$. However when could we be sure that the past innovations are more or less reasonably included into the stochastic trend? Is there any good proxy for random-walk's burn-in period? I've looked for R's burn-in default suggestion for $AR(p)$ part: `ceiling(6/log(minroots))`, meaning that in case of unit root we get infinity here and roughly tried to take close to unity root 1.0001 (actually it is the same like taking 60000 at once). So any reasonable suggestions or rule of thumbs you do use in practice?


title: Worst classifier

text: What is the worst classier that learns badly in practical problems? Edit: Especially bad on test data.. Thanks


title: What is the difference between 2x2 factorial design experiment and a 2-way ANOVA?

text: What is/are the difference(s) between 2x2 factorial design and 2-way ANOVA? Do they have the same assumptions?


title: How to analyze an unbalanced within-subject ANOVA design?

text: I have data from an experiment involving 4 groups of subjects, 2 possible interventions first and 3 possible intervention at a 2nd point, repeated data measurements from each subjects multiple time points with a decreasing number of measurements at each progressive time points after 2nd intervention (this was planned). My problem is that I am not sure how to analyze unbalanced data. More specifically, my experiment involves the following groups & interventions:   1. Control group (no 1st treatment, no 2nd treatment)     2. Experimental group 1 (1st treatment, no 2nd treatment)     3. Experimental group 2 (1st treatment, type A 2nd treatment)     4. Experimental group 3 (1st treatment, type B 2nd treatment)   The 1st treatment is the same for all subjects that received it. All the groups started out with an equal number of subjects (24 subjects) in each group (for a grand total of 96 subjects in the study). Samples were taken from the subjects at the following different time points to measure one dependent variable:   * `t(0)` \- Baseline AKA prior to time of the 1st treatment (or equivalent time) for ALL subjects (96 total subjects, so `96 total data points`)   * `t(1)` \- x number of weeks after the 1st treatment & prior to the 2nd treatment (or equivalent time) for ALL subjects (96 total subjects, so `96 total data points`)   * `t(2)` \- x number of weeks after the 2nd treatment (or equivalent time) for ALL subjects (96 total subjects, so `96 total data points`). 1/4 of the subjects from each group (6 subjects from each group, 24 subjects total) have the treatment site removed for analysis.   * `t(3)` \- y number of weeks after the 2nd treatment (or equivalent time) for remaining 3/4 of subjects (72 total subjects, so `72 total data points`). Additional 1/4 of the subjects from each group (6 subjects from each group, 24 subjects total) have the treatment site removed for analysis.   * `t(4)` \- z number of weeks after the 2nd treatment (or equivalent time) for remaining 2/4 of subjects (48 total subjects, so `48 total data points`). Additional 1/4 of the original number of subjects from each group (6 subjects from each group, 24 subjects total) .   * `t(5)` \- zz number of weeks after the 2nd treatment (or equivalent time) for remaining 1/4 of subjects (24 total subjects, so `24 total data points`). All the remaining subjects from each group have the treatment site removed for analysis. (I hope all that makes sense.) I am trying to analyze whether the 1st treatment has the effect we expect on the subjects, and whether the 2nd treatment has any further effects. Unfortunately, we _need_ to remove the treatment site after each timepoint `t(2)`-`t(5)` in order to obtain the microscopic data see what, if anything, is changing. The treatment area is quite small, so taking a portion of it is really not an option. Since I have 1 dependent variable to measure, 3 treatment variables, and don't have any random variables in my study, I am assuming I should use repeated measures fixed effect ANOVA. However, given I have progressively fewer data points (but equal among groups) for the last 3 time points (AKA starting with 96 @ `t(2)` \--> 72 @ `t(3)` \--> 48 @ `t(4)` \--> 24 @ `t(5)`), I am unsure how I should even start doing the analysis or what program would work for me. I have GraphPad Prism, but it only lets me go up to Two-way ANOVA, and I am unsure if I need to use 3-way ANOVA or not, or even if I need to use a different method. Suggestions on how to analyze the data and what program, if any, would be very helpful. Any other insights or suggestions are also welcomed.


title: Does Box-Cox parameter estimation count towards parameters for AIC?

text: Suppose I have a regression model with e.g. 2 parameters $y = ax + b$ But the data are non-normal so before regressing I transform both sides with Box-Cox estimation. Thus I get two Box-Cox parameters as well, $\lambda_x$ and $\lambda_y$. Now I want to calculate AIC for this model. How many parameters are there? My instinct would be that $\lambda_y$ counts as a parameter but not $\lambda_x$, because if the model is applied to forecasting $y$ from $x$, $\lambda_x$ can be estimated any time from the available $x$ but we have to remember which $\lambda_y$ to use as we don't know the $y$ we are trying to predict.


title: sum of discrete random variables - am I thinking right?

text: I have thirty discrete random variables for a risk management application. Each of the random variables may have a value $0$ (with probability $p_1$) or $X_i>0$ (with probability $p_2$). All thirty $X_i$ may be different. I am interested in the sum $S=\sum_{i=1}^{30} X_i$ of all thirty random variables. Am I right regarding the following thoughts: 1) an exact analytic solution for the probability distribution of S would demand the computation of $2^{30}$ different value-probability-pairs 2) the expected value of S can easily be calculated as the sum of 30 expected values 3) the standard deviation of S can only be calculated easily, if all 30 have a mutual covariance of 0 4) for a risk management application, which must consider values as well as probabilities, we must apply a Monte-Carlo-approch because of 1)-3) as there is no quick algorithm for an analytic solution 5) all would be even more clearly demanding the Monte-Carlo-Approach as in 4), if we have three differing instead of two differing values for the 30 random variables Before posting the question, I checked http://math.stackexchange.com/questions/295363/sum-of-two-random-variables which only speaks about identical discrete distributions. Othe questions I saw on cross-validated seem to address special discrete distributions like Bernoulli or continuous distributions like the normal Distribution.


title: multiple time series classification using randomForest R

text: I have daily time series of all constituent stocks/members of the S&P500, over a 5 yr horizon, and wish to classify as to whether a stock will report an "earnings revision" via a binary outcome (0=no earnings revision, 1 = earnings revision). I have been using a logistic regression, but would like to compare this against a machine learning method. Is the randomForest package in R capable of handling multiple time series for a classification problem? Any advice on getting me started would be very much appreciated. Thanks


title: Are there any techniques that quantify the importance/signification of individual attribute values of a particular data point?

text: Are there any techniques that quantize the importance of individual attribute values in a particular data point, in terms of the attribute's overall importance/signification/contribution to the uniqueness of the data point (uniqueness with respect to the whole data set). I'm not entirely sure if I've put it perfectly. Forgive if I'm not perfect in my presentation. I've looked up a lot on the web, particularly in Google Scholar, but couldn't find anything relevant while I'm running out of keywords I can think of. Please let me know if anything similar to what I've written exists.


title: Transforming time series of different time horizon to stationary

text: I have a list of monthly time series data with different time periods and different order of integration. I want to transform them all to stationary and a same time period. I noticed that the order of integration of a time series could be different with different number of data points. For example:   * TS$_1$ is I(1) with 67 obs.    * TS$_2$ is I(1) with 108 obs., I(2) with 67 obs., I(1) with 66 obs., and I(2) with 65 obs.   * TS$_3$ is I(1) with 108 obs., I(0) with 67 obs., and I(1) with 66 obs., etc. Suppose my data start at different times, but end at the same date. What is the best practice to transform them to stationary? Thank you!


title: What is the meaning of $A(B^2)$ interaction in a 3 level factorial with 2 factors $A$ and $B$, when the factors are qualitative?

text: I understand that when the factor is quantitative the $AB$ interaction can have a linear component and a quadratic component and so $AB^2$ makes sense there. But when we have qualitative elements how does this translate? _The Design and Analysis of Experiments_ by Montgomery does say that components $AB$ and $AB^2$ have no actual meaning. But why would we use them then? Can anyone provide me an intuitive answer to this or point to where I could find a better interpretation?


title: How to classify Gaussian from bimodal mixture of 2 Gaussians?

text: I have two types of samples:   * some look like bimodal Gaussian mixture (Type A),    * and some look like a Gaussian (Type B). How can I programmatically classify/label them? Language can be R, Python, Matlab, or anything else appropriate. ![enter image description here](http://i.stack.imgur.com/2zEVR.png) The data presented in the charts are Red values from jpeg images. In R code:               # I read an image     I = readJPEG(file, native = FALSE)     # so I is now a matrix of 3 dimensions (rows, columns, 3 for Red/Green/Blue)          # I extract a vertical line from the image, and only the Red part     image_extract <- I[150:260, 194, 1]          # After reading several images, I plot the 3 images image_extract for each type (A,B)     plot(image_extract_1)     lines(image_extract_2)     lines(image_extract_3)      For Type A I plotted, 3 image extracts on the same chart. Same for Type B. I hope it clarifies.


title: Neural net: variable distribution, validation and number of variables - unexpected results

text: I have been experimenting with ANNs on one of my datasets, they seem to have the potential to be quite effective in explaining my Y variation. Something i am finding is that they very much benefit from some prior modelling with other methods to reduce x variable quantity (i am finding i cannot use PCA as i need my X variables maintained and PCA prevents me doing this due to its factoring in this particular case not resembling anything in the real world). What appears to happen if i use all my variables is that some variables will be come out as powerful influences of my Y variable despite this absolutely not being the case in reality. The thing that bugs me about this is often the training and validation indicate the model is OK despite this, this may be due to not enough data points. The variables this occurs with are usually not very well distributed and typically i would remove these due to this. I am mainly asking as i was wondering why this happens? And also is it ok to prior model with other methods to remedy this? Also, how should i deal with validation if i a using two modelling types that relate to one another? Should the validation set be re generated each time? I have not used ANNs before (apart from out of intrigue) and do not fully understand exactly how they work in regards to re-correcting weightings through propagation


title: multiple imputation for a longitudinal study

text: i have an experiment wherein respondents were tested in two time points. however, respondents were tested at t1 and t2 OR t1 and t3 OR t1 and t4. Hence, data is missing at t2,t3,and t4 for 3/4 of respondents. i meet the missing at random assumption, but is there anything else I need to consider before i start up the imputation machine? thanks!


title: estimate conditional time-to-event for recurring event history model

text: I have a sample of individuals for whom I have measured the time at which they display a particular behavior. This behavior is recurrent (it is not fatal). I have modeled these events with time-dependent covariates, using piecewise constant hazard rates (based on an exponential distribution). This gave me parameter estimates for the effect of the different covariates; these parameters I assume to be constant over time. Now, knowing that an individual has values for covariates X1(t) and X2(t) at time t, I can calculate that person's fitted hazard rate at t (as exp(labda0 + labda1*X1(t) + labda2*X2(t)). The t's of interest are the points in time at which the behavior has been displayed. My aim is to calculate the time-to-event for a given individual as a function of the value for these covariates. So, given values for X1(t) and X2(t) and the parameter estimates for X1 and X2, I would like to estimate how long it will take take (in minutes) until the individual engages in the behavior next. How do I estimate that amount of time? It seems to me like a straightforward thing to calculate, but how?


title: Dawid Skene model for BUGS

text: I'm a beginner in BUGS. I'm trying to code the Dawid Skene model in BUGS. The model is as follows: ![](http://i.stack.imgur.com/hQRti.png) Currently, I'm adapting the code from Stan, here. Here is the code:               model {         for (i in 1:I) {          z[i] ~ dcat(pi[])        }        for (n in 1:M) {          y[n] ~ dcat(theta[jj[n],z[ii[n]]])        }     }       The data are               list(I=4, M=12, ii=c(1,1,1,2,2,2,3,3,3,4,4,4), jj=c(1,2,3,1,2,3,1,2,3,1,2,3),y=c(1,0,1,1,1,0,1,1,1,0,0,1))      I cannot compile the code since I keep getting the error: "variable pi is not declared in model or data". Do you have any idea how I can model the Dawid Skene model in BUGS? Edit: the context of this model is as follows: there are J annotators who need to label N items into K categories. There are M annotations ($M=N \times J$). The label of annotator $i$ for item $j$ is denoted as $r_{i,j}$. I need to find the actual correct label of each item $z_i$ and the competency of each annotator $\Theta$. Edit 2: $\Theta$ is a $J \times K \times K$ matrix where each element in the matrix is in $[0,1]$. In addition, if $i,j$ are fixed, $\sum_{k=1:K}\Theta[i,j,k] = 1$. $\rho$ is an array of $K$ elements. Each element in the array is in $[0,1]$ and the sum of the elements is 1. Edit 3: $r[i,j]$ is a categorical distribution taking $\Theta[i,j,:]$ as its parameter and returns an integer in $[1,K]$. $z[i]$ is also a categorical distribution which takes $\rho$ as a parameter and returns an integer in $[1,K]$. The question is that I don't know how to model $\Theta$ in BUGS. Since it is also a parameter that I want to infer, I suppose I need to model $\Theta$ as a matrix in the model (which I don't know how).


title: Finding categorywise significance in regression

text: I have a survey dataset which has "Overall quality of experience" as the dependent variable. I further have 4 sub-categories related to Food, Menu, Facility and Team. Each sub-category has some questions(5,5,3,3 respectively). Note that none of the sub-categories has an overall question such as "How is the overall quality of food" or "The overall quality of menu". After running an Ordinal Logistic regression, I am able to arrive at the individual significance for each of the 16 independent variables. However I want to calculate the overall individual significance for 'Food, Menu, Facility and Team'. Could you guys help me out here? Thanks!


title: Estimate linear regression using items randomly selected from an item pool

text: I am asking this question against the background of a linear regression with single predicted variable $Y$ and multiple predictors $X$. $X$ comes from a survey using an "item pool" which suggests that not all items are presented to each respondents, but only a subset of items is. This subset is randomly selected (but in the current case $Y$ is always presented to all respondents). The objective is to estimate a linear regression model on all $X$ irrepesctive of the missing information due to randomized selection of items. I have the following three questions: 1) Given randomization it seems to me that the data are missing completely at random (MCAR). Thus, estimates of sufficient statistics (correlation matrix and means) are unbiased. Hence, I infer that regression coeffcients should be unbiased. Is that correct? 2) Standard software will listwise delete cases resulting in an empty set of observations. Should I estimate coeffcieints using the sufficient statistics only? What would be suitable approaches in R (e.g., available functions)? 3) To add complexity to the problem, assume items vary in (known) probability to appear in a survey. Thus sample size varies per pair of items as a function of $np_1p_2$, where these are sample size and probabilities to appear in a survey for any two items. How should the varying sample sizes be included in estimating standard errors of estimate, and again how can available functions in R deal with this?


title: Cross-Validation probability

text: I'm fitting a bulge-disc sersic model to galaxy light profiles and need to find a robust method of finding the errors on the parameters. I'm looking at cross-validation implemented in this paper here (page 9, section 4.2). It says I need to calculate the likelihood that the point has been left out. How do I compute this? I'm not a statistician and my online research has only confused me! I've never looked at likelihoods before!


title: Sum of random vectors with fixed amplitudes

text: Is there a simple way to evaluate the pdf of the amplitude of a sum of vectors with fixed amplitudes and random phases? Explicitly, let $Ae^{i\phi}=\sum_{n=1}^NA_ne^{i\phi_n}$, where $N$ is a fixed positive integer, the $A_n$ are fixed positive real numbers, and the $\phi_n$ are independently uniformly distributed over $[0,2\pi)$. Then $A$ is the non- negative real magnitude of the sum, and $\phi$ is a phase that I don't care about. Is there a simple expression for the pdf of $A$? (Or equivalently, a simple expression for the pdf of $A^2$?) Failing that, is there a reasonable approximation that works for both large and small $N$? I'm working from this survey of research on related questions:   * Abdi et al. 2000, "On the PDF of the sum of random vectors", IEEE Trans. Commun. **48** (1) 7-12. DOI: 10.1109/26.818866, free preprint: http://web.njit.edu/~abdi/srv.pdf I've been wading through the references given on the first page for item (3) and (c) "various analytic approximations", but I haven't found anything satisfactory. It appears to be a hard problem in general. The good news is that I'm interested in fixed amplitudes $A_n$, which should simplify the problem. The bad news is that I'm interested in arbitrary values of the $A_n$. There could be 100 vectors of similar magnitudes, or there could be 1 dominant vector and 5 much smaller vectors, etc. I mention this because some references treat only the case where all of the $A_n$ are equal. I want to evaluate the pdf millions of times per second, which is one reason why I'm shying away from the integrals in part II of the survey. If the integrals undergo a massive simplification for fixed $A_n$, it's not obvious to me, so please let me know if that's the case! I'm hoping to use the answer to this question in a Bayesian approach to polyphonic pitch tracking. It would be useful to evaluate the likelihood of an observed filter response arising from a proposed set of tones with specific filtered amplitudes and unknown phases. If there's no great solution, I could fall back to an ugly hack: If $N=2$ then compute the exact solution; if $N>2$ then approximate the distribution of $A$ with a Rayleigh distribution, truncated to the correct support. So, I guess I'm looking for a solution that's a little better than that. For example, it would be great to put peaks in the correct places when $N=3$.


title: Using principal component analysis (PCA) for feature selection

text: I'm new to feature selection and I was wondering how you would use PCA to perform feature selection. Does PCA compute a relative score for each input variable that you can use to filter out noninformative input variables? Basically, I want to be able to order the original features in the data by variance or amount of information contained.


title: Dimension reduction using Principal Component Analysis

text: I have basic question regarding PCA. I am given a task of reducing independent variables in a dataset of about 250,000 enteries and 400 independent variables because some of them are highly correlated. So I want to retain only the most useful ones. I want to perform PCA analysis on it and I have read up but I still don;t see how it results in dimension reduction. I perform the analysis and select the Principal Components based on several criterias. Now what? How do I chose the most relevant of those 400 variables and delete the rest? Thank you.


title: Bayesian parameter estimation homework problem

text: For BPE approach,for p($\theta$)~U(0,1) show that p(x|D)= $1/n+2[(x+ \sum_{k=1}^n(x_k)!)$ $ (n+1-x-$($\sum_{k=1}^n(x_k)$!]$/$[$(\sum_{k=1}^n(x_k))$! $(n-\sum_{k=1}^n(x_k))$!] I have tried using the following equation p(x|D)= $\int p(x|\mu)p(\mu|D)d\mu$ I have substituted for the 2 terms but unable to reach the required solution


title: Different degrees of freedom when using the same mixed effect model (SPSS)

text: I have a question regarding mixed effect models. Is it possible that the degrees of freedom differ when I use the same model with different dependent variables? Or does that indicate that something went wrong? I analyzed my data using the MIXED procedure in SPSS. There are two random factors (participants and stimuli) and one fixed (continuous) factor (PD) in my design. The number of degrees of freedom is 1201,436 for my first dependent variable "perf". When I use the same model and only choose a different dependent variable ("wr") the number of degrees of freedom is 1194,166. The number of data sets and the model dimensions are identical for both dependent variables. The information criteria differ (e.g., AIC for "perf" = 3975 and AIC for "wr" = 9942). The tests of fixed effects are significant for both dependent variables. Does it make sense that the degrees of freedom differ?


title: Can someone sort me out regarding the calculation of AUC?

text: I am having some trouble with two different implementations of a classification problem giving different results. Me and my college who did the other implementation has narrowed the problem down to the way we calculate the area under the receiver operating characteristic curve (AUC). One solution is derived from a formula appearing at least at one location: [1] $$AUC_1 = \frac{1}{mn}\sum_{i=1}^{m}\sum_{j=1}^{n}\mathbf{1}_{p_i<p_j}$$ I have ported the implementation into `R` and compared it with the result $AUC_2$ from the R-package `pROC`:               auc1 <- function(p) {       values <- p[,1]       positives <- values[ p[,2]==1 ]       negatives <- values[ p[,2]==0 ]       count <- 0       for ( p in positives ) {         for ( n in negatives ) {           if ( p>n ) {             count <- count + 1           }         }       }       return(count/(length(positives) * length(negatives)))     }          auc2 <- function(p) {       library(pROC)       c <- roc(p[,2], p[,1], print.auc=TRUE, ci=F, of="auc")       return (auc(c))     }          predicted <- c(0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1)     real      <- c(0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1)     p <- cbind(predicted, real)          auc1(p)     auc2(p)      Can someone shed some light on why I get: $$AUC_1 = 0.9090909 \mbox{ and } AUC_2 = 0.9545\mbox{?}$$


title: How do I interpret outcomes of principal component analysis?

text: At the moment I am learning to use principal component analysis (PCA). The whole idea of the method is quite obvious to me, and on this I have no questions. However, I do need some explanation on how to interpret the output of my computations. I have seen posts, Interpreting PCA scores and What are principal component scores?, which helped me a lot, but still I am unable to "get something out" of my data. I have created a dataset consisting of variables $x$, $y$, $z$. Variable $z$ depends very much on $y$, while $y$ and $z$ almost do not depend on $x$ variable. I have performed PCA computations, from which I've gotten eigenvalues and eigenvectors (sorted in the direction of decreasing eigenvalues). My question is: How can I conclude which variables are meaningful and which are meaningless?   How can I say which variables are correlated with each other? What else could I possibly squeeze out of the data, and how? EDIT: I am asking what do I do with the modified data? I have new set of data which has less dimensions. What now? How does it help me to understand the data?


title: Standard deviation for adjusted effect size estimates

text: I am trying to calculate the effect size for a trait adjusted for my set of nuisance variables. My models looks like this:               Trait ~ Dx + Age + Sex      Where "Trait" is a continuous variable and "Dx" is a binary group membership variable. I'm using this basic formula to get an effect size estimate:               (mean(Trait[which(Dx==0)])-mean(Trait[which(Dx==1)]))/SDpooled      but I am mostly interested in the effect size after correcting for `Age` and `Sex`. I calculated the marginal means (i.e. corrected for Age and Sex) for the `Trait` split by `Dx` group as shown here:  on page 15, section 7.5 (PDF) But I want to know which standard deviation I should use. Is it required to also adjust the pooled standard deviation by the covariates? Is there some better way to get effect size estimates adjusted by a set of nuisance variables?


title: How to interpret a regression coefficient for the reciprocal of an independent variable?

text: Does anyone know how to interpret a coefficient when the variable in the model is the reciprocal of the original variable? I have an inverse equation, where $\text{time} = \beta_0 + \beta_1(1/\text{horsepower})$, where $\text{time}$ is how long it takes to accelerate a car. I need to hypothesize the sign of the variable $\text{horsepower}$. But to do that, I need to understand what it means.


title: When sampling distribution and posterior distribution under uniform prior belong to the same family?

text: This is similar to conjugate prior, except that here I consider if it is possible that a sampling distribution and the posterior distribution under uniform prior belong to the same family of distributions? In other words, when viewing $p(x|\theta)$ as the likelihood function of $\theta$ given $x$, and then normalize it to be a distribution for $\theta$, are there examples where $p(x|\theta)$ and the distribution for $\theta$ belong to the same family? More generally, are there examples for priors other than uniform and sampling distributions, where the posteriors fall into the same family of distributions as the sampling distributions? Thanks!


title: Model selection with Firth logistic regression

text: In a small data set ($n\sim100$ ) that I am working with, several variables give me perfect prediction/separation. I thus use Firth logistic regression to deal with the issue. If I select the best model by AIC or BIC, should I include the Firth penalty term in the likelihood when computing these information criteria?


title: Determine the optimum learning rate for gradient descent in linear regression

text: How can one determine the optimum learning rate for gradient descent? I'm thinking that I could automatically adjust it if the cost function returns a greater value than in the previous iteration (the algorithm will not converge), but I'm not really sure what new value should it take.


title: Statistic test on percentage correct classified by emotion recognition

text: For a potential emotion recognition bachelor-project I was wondering what statistical test I have to perform when I get my results to test whether it's significant. I will be testing which combination of feature extraction and machine learning algorithm will give me the best percentage of correct classified. The results will exist out of Combination A gives ...% classified, Combination B ...%, Combination C ...% and so on. Which statistical test should I use to test whether Combination ... is significantly better than the others and why? For example: 6 emotions have to be recognized in a database with 100 faces for each emotion (600 total). Every machine learning algorithm will use 2/3 for training and 1/3 for testing. Which face per emotion is in the training set and which one is in the test set is randomly selected every epoch for 100 epochs. The end result is for example: Combination A classified 93.4% correct, Combination B 91.2%, Combination C 86.3% and so on. Which statistical test should I use to test whether Combination A is significantly better than Combination B (and C) and why? Also is it dichotomous? As the probability of successful selecting the right emotion is 16.67%.


title: Excluding factors

text: I want to calculate the influence of different factors on a list of y. The problem is, that some of the y are also include in the factor matrix, so cv.glmnet would compute a coefficient around 1 for this factor, but I want that it is zero, so only the other factors are considered for the coefficient matrix. So my question is, are there any options to exclude the factor in this case or to change the coefficient to zero during the calculating? My easy solution would be to set the values of the factor to zero, so that it can't have an influence. Script:               for ( i in 1:dim(allsamples)[1])     {        #allsamples= x genes, Y samples        gene=as.data.frame(allsamples[x,])             #for one gene --> Y observations        gene=as.numeric(gene)        gene=gene[-1]             #expression of the influence factors        tfs=t(data.matrix(tfgenes))        colnames(tfs)=tfgenes[,1]        tfs=tfs[-1,]             # set expression to zero if gene=influence factor?       if(length(intersect(gene,colnames(tfs)))==1)       {           tfs[,intersect(gene,colnames(tfs))]=rep(0,125)  #row or column??       }else       {           tfs=tfs       }      z=NULL           #repeat regression            for(i in 1:100)           {              cv.fit=cv.glmnet(data.matrix(tfs),gene,family="gaussian",nfolds=10,                               alpha=1,type.measure="mse")              lambda=cv.fit$lambda.min              z=rbind(z,lambda)          }           # all lambda.min      mean_lambda_min=mean(z)           #regression with the mean lambda minimum      fit=glmnet(data.matrix(tfs),gene,family="gaussian",lambda=mean_lambda_min,alpha=1)           # matrix of coefficients for each gene        matrix=rbind(matrix,t(as.matrix(coef(fit))))     }


title: Using MANOVA for classification without separating training and test sets

text: In this study: Rosenblum, Sara, et al. "Handwriting as an objective tool for Parkinson’s disease diagnosis." Journal of neurology 260.9 (2013): 2357-2361 The researchers attempt to classify Parkinson's disease (PD) with a control group of 20 participants and a PD patients group of 20 participants, while the participants are writing on a digitized pen and paper. The researchers found that a MANOVA analysis predicts which is which with 95% accuracy. I'm pretty new to MANOVA, but it seems to me that there is a need to separate the groups to training set and validation sets. Am I correct? In other words, is MANOVA prone to overfit?


title: What is a good way to implement inverted LRT confidence interval?

text: One standard method to derive the confidence interval of a parameter is the LRT method, in which we try to find the acceptance region of the LRT test with size $1-\alpha$, then "inverting" it we get a confidence interval of confidence $1-\alpha$. A proof of it can be found in Berger & Casella, page 421. However, in practice implementing this method is not easy even for simple distributions like the Poisson distribution, the Binomial distribution, etc. The strategy is: For each $\theta_0\in \Theta_0$, let $A(\theta_0)$ be the acceptance region of a level $\alpha$ test of $H_0:\theta=\theta_0$. For each $x\in \mathcal{X}$, define a set $C(x)$ in the parameter space by $$ C(x)=\\{\theta_0:x\in A(\theta_0)\\} $$ Since in practice we do not really know what $\theta_0$ is, and has to reply on various estimates like MME or MLE, the above strategy is computationally intensive. To implement it in R, we would need to select an open interval $O$ contain $\theta$, randomly choose $N$ points in $O$, then for each point $\theta_{i}$ find the acceptance region and whether $x$ is in it (usually numerically by evaluating over $M$ points in $\mathcal{X}$). Finally we choose the smallest and largest $\theta_i$ as the confidence interval. The complexity lies in the LRT statistic. For each $\theta_i$, to find the acceptance region we first need to compute $c_{\alpha,i}$. Then "invert" the LRT test to find the acceptance interval $\lambda(x)>c_{\alpha,i}$. And this is not trivial as we have to compute a valid $c_{\alpha,i}$ that determine a rejection region with level $\alpha$. In practice, even for binomial distribution this is not easy (though it can be done). And for a more complicated distribution it might be a very complicated process. This, coupled with the fact that we need to evaluate over $N$ points, leads to the computational complexity exceeding $O(NM)$ or maybe even $O(NM^2)$. Further, taking approximations at every step in the numerical computation would make the actual result less accurate. So how do we actually do this in practice? Has there been any effort to reduce this computational complexity? I ask here as I could not find an easy answer from elsewhere. I assume others must have asked the same question earlier.


title: Continuous Time Markov Chain Transition Rates

text: > A hospital has two physicians on call, Dr Dawson and Dr Baick. Dr Dawson is > available to answer patients' calls for time periods that are exponentially > distributed with mean 2 hours. Between those periods, he takes breaks, each > of which being an exponential amount of time with 30 mins. Dr Baick works > independently from Dr Dawson, but with similar work patterns. The time > periods she is available to take patients' calls and the times she is on > break are exponential random variables with means 90 and 40 mins, > respectively. In the long-run, what is the proportion of time in which > neither of the two doctors is available to take patients' calls? The problem I am having with this question is not finding the steady state probabilities, but rather, how to set up the transition rate matrix for this question. My approach is as follows: Let $X(t) = 0$ if neither Dr Dawson nor Dr Baick is available to answer patients' calls. Let $X(t) = 2$ if both are available to take calls. Let $X(t) = d$ if Dr Dawson is available to take calls and Dr Baick is not. Let $X(t) = b$ if Dr Baick is available to take calls but Dr Dawson is not. Then $\\{X(t), t\ge 0\\}$ is a continuous-time Markov Chain with state space $S=\\{0,2,d,b\\}$. Let $q_{ij}$ denote the transition rate from state $i$ to state $j$, for $i,j \in S$. Then it can be shown that $Q$ (the matrix of all $q_{ij}$'s) is given by: ![enter image description here](http://i.stack.imgur.com/hwrjg.jpg) where the rows are labelled in the order: $0, 2, d, b$ and same with columns. What I'm confused is that how is $q_{20} = 0$? Isn't it possible to transition from the state where both doctors being available to answer calls to no doctors are available?


title: How to calculate summary statistics for nested factors?

text: I have a data set where I have two levels of within variables. I have a set of `Scenarios` and each scenario has a set of `Trials`, and all scenarios are played by the same people. As I understood from here my `Trials` are nested in the `Scenario`s. In this case how can I compute the means and standard deviations? And to take a step further if I have multiple sets of people -- in independent categories, like in the example from the reference, how can I specify the ANOVA model which describes my data? (Say the sets of people play in different `conditions`, but the same `Scenarios` with the same `Trials`).


title: Variance of Z for Z = X + Y, when X and Y correlated

text: So I'm trying to show that ${\rm Var}(Z) \le 2({\rm Var}(X)+{\rm Var}(Y))$ for $Z = X + Y$. This seems to be pretty easy to show given that $X$ and $Y$ are uncorrelated. But I'm running into trouble at this step: $$ {\rm Var}(Z) = {\rm Var}(X) + {\rm Var}(Y) + 2E[XY] - 2E[X]E[Y] $$ Normally you could say, $X$, $Y$ uncorrelated $\rightarrow E[XY] = E[X]E[Y]$, but when you _cannot_ do this, I'm lost. Any tips?


title: Quadratic (programming) Optimization : Multiply by scalar

text: I have two - likely simple - questions that are bothering me, both related to quadratic programming: **1).** There are two "standard" forms of the objective function I have found, differing by multiplication of negative 1. In the R package `quadprog`, the objective function to be minimized is given as $-d^{T}b + \frac{1}{2}b^{T}Db$ and in Matlab the objective is given as $d^{T}b + \frac{1}{2}b^{T}Db$. How can these be the same? It seems that one has been multiplied through by a negative 1 (which as I understand it would change from a min problem to a max problem. **2).** Related to the first question, in the case of using `quadprog` for minimizing least squares, in order to get the objective function to match the standard form, it is necessary to multiply the objective by a positive 2. Does multiplication by a positive number not change the solution? **EDIT:** I had the wrong sign for the Matlab objective function.


title: Do I have to learn SAS if I want to go into industry?

text: I am currently a PhD student in Biostatistics (finishing my 2nd year). My original goal when entering graduate school was to go into academia. Recently I have been debating this choice and wondering if perhaps industry will be a better fit. I come from a strong programming background. I love using R because it is a lot more similar to the programming languages that I grew up with, and is even able to interface with them. I abhor SAS. Whenever I am forced to use it, I feel like I am wearing 2 left shoes that are 3 sizes too small. I know that SAS is predominant in industry however. Do you guys think I will just have to suck it up and learn SAS proficiently, or will I be able to get by with R?


title: Poisson regression with large data: is it wrong to change the unit of measurement?

text: Due to the factorial in a poisson distribution, it becomes unpractical to estimate poisson models (for example, using maximum likelihood) when the observations are large. So, for example, if I am trying to estimate a model to explain the number of suicides in a given year (only annual data are available), and say, there are thousands of suicides every year, is it wrong to express suicides in hundreds, so that 2998 would be 29.98 ~= 30? In other words, is it wrong to change the unit of measurement to make the data manageable?


title: Is it valid to judge the superiority of a classification technique using only artificial data?

text: ... where the significance of the results has been checked with statistical tests. If no, then why? **side note:** The data is generated based on specified prior probabilities.


title: Should I run separate regressions for every community, or can community simply be a controlling variable in an aggregated model?

text: I am running an OLS model with a continuous asset index variable as the DV. My data is aggregated from three similar communities in close geographic proximity to one another. Despite this, I thought it important to use community as a controlling variable. As it turns out, community is significant at the 1% level (t-score of -4.52). Community is a nominal/categorical variable coded as 1,2,3 for 1 of 3 different communities. My question is if this high degree of significance means I should be doing regressions on the communities individually rather than as an aggregation. Otherwise, is using community as a controlling variable essentially doing that?


title: How to deal with outliers?

text: **Background:**   There are often some values among sampling set that appear not closely compatible with the rest. They would be called as extreme values or simply outliers. Dealing with outliers has been always a matter of challenge. There are some approaches to solve the problem of the existence of outliers:   * moving them to a separated set   * replacing them with nearest values from non-outlier set   * ... **Question:**   What is the most recommended method(s) to deal with outliers? (with details and an example)


title: Reason for taking the mean of experimental observations

text: Whenever we do some experiment in Physics, e.g. measure Planck's constant, we take the mean of several experimental observations. This seems to be the most obvious way to report a result. Is there a reasoning based on MLE or other statistics concepts that can be used to justify that the mean is the best way to measure the value over several observations?


title: How large should a sample be for a given estimation technique and parameters?

text: Is there a rule-of thumb or even any way at all to tell how large a sample should be in order to estimate a model with a given number of parameters? So, for example, if I want to estimate a least-squares regression with 5 parameters, how large should the sample be? Does it matter what estimation technique you are using (e.g. maximum likelihood, least squares, GMM), or how many or what tests you are going to perform? Should the sample variability be taken into account when making the decision?


title: Spatial statistics models -- CAR vs SAR

text: When would one prefer to use a Conditional Autoregressive model over a Simultaneous Autoregressive model when modelling autocorrelated geo-referenced areal data?


title: How to deal with the effect of the order of observations in a non hierarchical cluster analysis?

text: When a non-hierarchical cluster analysis is carried out, the order of observations in the data file determine the clustering results, especially if the data set is small (i.e, 5000 observations). To deal with this problem I usually performed a random reorder of data observations. My problem is that if I replicate the analysis n times, the results obtained are different and sometimes these differences are great. How can I deal with this problem? Maybe I could run the analysis several times and after consider that one observation belong to the group in which more times was assigned. Has someone a better approach to this problem? Manuel Ramon


title: Cook's distance and $R^2$

text: I am currently running a linear regression and calculating its $R^2$ After that. I calculate the Cook's distance of all points and throw away from the analysis all of the points with a distance higher than $d_i >\frac{4}{\text{No. observations}}$. To my surprise the $R^2$ is worse. How is this possible?


title: Confusion about hidden Markov models definition?

text: In Wiki the following diagram has been sketched and the following text has been written: ![enter image description here](http://i.stack.imgur.com/Jr13M.png) > From the diagram, it is clear that the conditional probability distribution > of the hidden variable x(t) at time t, given the values of the hidden > variable x at all times, depends only on the value of the hidden variable > x(t − 1): the values at time t − 2 and before have no influence. Now the problem is I don't understand why this is true for $X_1$ as there is a loop there and as a result I cannot see why this follows: $$p(X_1|X_2,X_3)=p(X_1|X_0=\emptyset)=p(X_1)$$ But it actually should be: $$p(X_1|X_2,X_3)=p(X_1|X_2)$$


title: Weighted regression for categorical variables

text: I have been trying to use weighted regression on some data where there are three categorical variables with the               lm(y ~ A*B*C)       command in R. This was to try and deal with the problem that my residual plot does not have constant variance, and that for one of the factors the variance for two of the levels is much smaller than the other level and in the other factors, so I define weights, in a vector ww, for each observation, so that the corresponding observations do not contribute as much as the other observation, and use               lm(y ~ A*B*C, weights=ww)      My problem is that I get exactly the same residuals when I use weights than when I don't use weights. Is this supposed to happen? Is it because all the variables A, B, C are categorical? It sort of seems to make some sense because the regression coefficients (the beta's in the linear model) are just the means and the differences from the reference mean in the un-weighted model, and I think that will also be the case using weighted regression, but I am not sure. Can anyone help me with this? I am a bit lost. i.e. can someone tell me whether I should NOT be using weighted least squares to fix problems with heterogeneous variance of residuals? (I have also tried transforming the response variable, $\log(y)$, it helped but didn't fix the problem.


title: Aggregating Survey Information Over Time

text: Suppose you have a monthly survey that asks a variety of questions that you aggregate the answers and standard error. However, some respondents begin to provide answers to the monthly survey on a more frequent, but irregular basis. Note that it is some, but not all of the respondents, and their answers will come on different days, potentially with different overlap. What is the best to update the survey on a daily basis to provide aggregate survey information (or a time-series of the underlying question response)? Replacement? Moving Average? Kalman Filter? Other techniques? Alternately, consider a situation where you have multiple surveys asking similar enough questions. You obtain information about the surveys on different days (means, counts, standard errors), potentially more than one a day. What's the best way to aggregate the different survey information to one underlying series?


title: Logistic regression goodness of fit

text: I have fitted a logistic regression model showing the probability of flight given distance to a disturbance event. My sample size is 140 observations of which 45 were observed to fly. Distance to the disturbance is the only significant predictor among many explored (biologically relevant interactions as well). Here is the out from the R console               Call:     glm(formula = flight ~ approach_km, family = binomial)          Deviance Residuals:         Min       1Q   Median       3Q      Max       -1.1233  -0.9317  -0.7460   1.3538   1.9400            Coefficients:             Estimate Std. Error z value Pr(>|z|)       (Intercept)  -0.1233     0.3438  -0.359   0.7199       approach_km  -0.4530     0.2225  -2.036   0.0418 *               Null deviance: 175.82  on 139  degrees of freedom     Residual deviance: 171.30  on 138  degrees of freedom     AIC: 175.3          Number of Fisher Scoring iterations: 4      Here is a plot of the predicted probabilities of flight. ![binary logistic regression, flight~distance](http://i.stack.imgur.com/V60uB.jpg) As you can see there is no point along the distance continuum in which the probability of flight exceeds 0.5. Since this is the case the classification tables will never predict a flight (it predicts 67.9% correctly, ALL THE NO FLIGHT EVENTS!). McFadden's psuedo R2 = 0.026. I am happy with the fit of the model but I am under the belief that I need to include these types of goodness of fit statistics when publishing. I would rather show that within a certain range of the data, say 1 km, there is roughly a 0.40 probability of flight. If there were 20 observations associated with this estimate there would be 8 flight events and 12 non-flight events. Is there a way to show this? How can I write up these results to show my findings without referee's balking at my lack of fit tests? Any thoughts or references would be greatly appreciated. Thanks. UPDATE: In response to @Drew75. Here is the ROC curve. Does this mean that the model does not perform as well I has thought? ![ROC curve](http://i.stack.imgur.com/p8uqj.jpg)


title: Non-parametric 2-sample test for correlated data

text: I have two samples of data, S1 and S2, that both suffer from spatial autocorrelation. The samples are different sizes and are not paired (specifically, S1 is comprised of M subsamples of k spatially correlated data points, while S2 contains N of these subsamples). The data are mostly continuous, however, due to the quantity being measured, there are some values which are undefined but known to be high compared to the rest of the sample, and have been set to fixed artificially high values. For this reason, only median statistics can be used and whitening the data using the covariance matrix is not possible. I would like to check if S1 and S2 come from the same population, however I am unable to find something non-parametric that works with correlated data. I was wondering if anyone is aware of a test that could be applicable in this situation.


title: Repeated measures 3-way ANOVA: main effects

text: I have 6 anatomical landmarks on one bone (humerus) and on 2 points on another bone (ulna). I measured distances between every different combination of those points in the range of motion (every 10º from 10º to 90º) of those 2 bones (elbow). I end up in a $2\times6\times9$ repeated measures layout. Basically I have 6 combinations over flexion angle for 1 point on the ulna, and another 6 combinations over flexion angle for the other point on the ulna. And my question is if there is an overall difference between the first 6 combinations and the other 6 combinations. Is it possible to answer this question with a 3-way repeated measures ANOVA, taking the main effects? I already compared each combination with the post-hoc analysis, but it would also be nice to know if there is an overall effect. I would really appreciate your expertise on this topic. Thank you so much in advance.


title: When accuracy in the cross-validation process less,is reducing the features a good idea?

text: I am doing a project for classifying the presence of cars/bikes in an image.I have extracted the features from the images(data-set of cars and images not belonging to that of cars) and applied K-means clustering to get a uniform feature vector X for all images.Now upon applying Cross Validation( Folds=10) upon the data-set yields me an accuracy of ~61% when trained using SVM classifier( RBF Kernel).In this case,is reduction in the number of features a good option(before giving the feature obtained from image before giving it to clustering)? Thanks and regards


title: If logit and neural net are very close?

text: What do you conclude if your logit and neural net results are very close in terms of precision. Where the logit has each variable going in independently. Should you conclude that your additively separable model in the logit is probably correct as it is close to the neural net result, which can approximate any functional form? Thank


title: Decision tree for censored data

text: Does anyone know an implementation in R (or other) of a decision tree for censored outcomes? I would like to utilize a decision tree in order to discretize/bin continuous variables before a survival analysis in some sort of principled manner. I am left with only a traditional decision tree using a binary target (event/no event) disregarding the censored nature of the data as it stands.


title: Proof for “Least squares estimator is BLUE”

text: This question was raised before. But so far I haven't got any satisfactory answer yet. So I raise it again. Suppose we have the linear model $$ y = X \beta + v, \quad v \sim (0, \sigma^2 I), $$ where $X$ has full column rank. The least squares estimator $\hat{\beta}=X^\dagger y$ is claimed to be the best linear unbiased estimator of the true parameter vector $\beta$. I checked all the books and on-line materials I could find for the proof, but found all of them seem to have a derivation problem. To prove it, one assumes $\bar{\beta} = Cy$ is any linear unbiased estimator of $\beta$. Using the fact that $\bar{\beta}$ is an unbiased estimator, we can easily obtain $(CX-I)\beta = 0$. Then all the books and on-line materials just conclude that $CX-I$ must be 0. I don't understand this at all. If $\beta$ were arbitrary, then certainly we would have $CX-I = 0$. But $\beta$ here is a specific parameter vector. Why can one claim $CX-I=0$ ?


title: Odd problem with a histogram in R with a relative frequency axis

text: Using either truehist() from MASS or just the normal hist() function in R with the prob=TRUE option, I'm getting very strange values for the y-axis. I was under the impression that these values should all be below 1.00, as the relative frequency of any value should be below 1.00 and the area under the curve adds to that. Instead, I'm getting axes with ranges toward 1500, and step sizes in the hundreds. Does anyone know what's going on? The values aren't event consistent, so it doesn't seem like they've got any relative scaling to them. For reference, I'm using the following code: `hist(g1$Betweenness, main="", xlab="Betweenness", sub="Generation 1", prob=TRUE)` The data for one such plot: 0.009619951 0.009619951 0.006750843 0.006750843 0.006750843 0.006750843 0.014497435 0.006750843 0.006750843 0.006750843 0.006750843 0.006750843 0.006750843 0.006750843 0.006750843 0.006750843 0.006750843 0.006750843 0.006750843 0.006750843 0.006750843 0.006750843 0.006750843 0.006750843 0.008663582 0.008663582 0.006750843 0.012058693 0.012489059 0.024587132 0.084941213 0.01248905 0.012489059 Annoyingly, JMP handles this just fine, but I've come to prefer R's plotting style.


title: Help in validating my experiment design

text: I've devised an experiment plan for an experiment in the domain of cognitive science, but I'm not 100% sure, if my idea about statistical analysis is correct. I'm performing this experiment for my Master's Thesis. My promoter is more of a philosopher than a statistician, so I'm turning to you as people experienced in the field. ( **Question 1** :)Could you please take a look and tell me, if it sounds ok or not? Additionaly I have two other questions, if you don't mind me asking :) # The experiment There are two groups of stimuli:   1. Images      1. of type _A_     2. of type _B_   2. Sentences      1. of type _x_     2. of type _y_ A person is first given a stimulus (an image) of type _A_ or _B_ , and then a second stimulus (a sentence) of type either _x_ or _y_. The person then has to indicate whether the sentence is correct for the preceding image. So I have 4 possible scenarios:   1. _A_ -> _x_   2. _A_ -> _y_   3. _B_ -> _x_   4. _B_ -> _y_ I would like to know:   1. _if_ there are interesting differences between, say, (1) and (3).   2. _the direction_ of the differences between cases (for example, people gave more correct responses in case (1) than in case (3)).   3. _the magnitude_ of the differences between cases. # Statistics ## Intro Factors, or indedepndend variables, are _I_ for images, and _S_ for sentences. Factor _I_ has two levels ( _A_ and _B_ ), and also factor _S_ has two levels ( _x_ and _y_ ). The outcome / dependent variable is the number of correct responses. The design is a _I_ x _S_ (two-way) factorial design. Each subject is given the whole spectrum of stimuli, so the design is a within-subjects, repeated measures. ## Results After checking if the neccessary assumptions (like homogeneity of error variances) had been met... In order to know:   1. _if_ there are interesting differences between, say, (1) and (3), I will first need to check if there is interaction between _I_ and _S_. If so, then I'll perform a post-hoc test (Tuckey's) to find out where exactly are the differences.   2. _the direction_ of the differences between cases (for example, people gave more correct responses in case (1) than in case (3)), I will see which mean is bigger on the correct position in the post hoc test.   3. _the magnitude_ of the differences between cases. Effect size is $R^{2}$ or $\eta^{2}$. ( **Question 2** :)But how to see it for individual cases? The last thing I don't know is how to interpret, or correct the situation, if it so happens that, for example there will be difference (main effect) between _I_ (images). I would rather have no differences between images and sentences, only the interaction. Of course I cannot count on that. ( **Question 3** :)Is there any way of tackling such problems?


title: Statistical meaning of pearsonr() output in Python

text: I have calculated the Pearson correlation by using               pearsonr(var1, var2)      I understand that the first number is the Pearson correlation and the second number is the significance. I have a few questions:   * Above which value can we consider significant correlation?   * Is the R-square only `R**2`?   * How do I calculate the adjusted R square?


title: Likelihood ratio test versus AIC for model comparison

text: When performing model comparison, why does the likelihood ratio test require two nested models while this is not required when using the AIC?


title: Linear vs polynomial regression

text: If you have multiple dimensions to your data, where it is not possible to visualize them together, how to decide if your model should be linear or polynomial?


title: How to find lambda in R, if other parameters are known

text: I am solving some problems on Poisson distribution and I have come across a problem where in I need to check the rate (lambda) of a particular thing with a probability of 0.001 (1 in 1000) of p(4). Specifically, I need to figure out the probability of finding a 1 cookie out of 1000 that has less than 5 raisins. So I need say how many raisins should be present in dough before its mixed to ascertain the above probability. **How am I suppose to do this in R?** Here is the actual question at hand (I am interested in the **Q.7a** ), and here is the link with answers, please look at the answers. > **7.** A batch of cookie dough will be sliced up into 100 cookies and then > baked. 400 raisins have been included in the batch of dough, and the dough > has been thoroughly mixed so as to randomize the ingredients. What is the > chance that, despite these precautions, one or more cookies in the batch > will contain no raisins at all? > > **7a. How many raisins should be put in the batch of dough to be 99% sure > that no cookie comes out with no raisins in it?**


title: What is the lasso in regression analysis?

text: I'm looking for a non-technical definition of the lasso and what it is used for.


title: ARIMA double seasonality with dummy in R error xreg

text: I'm playing with hourly univariate data and try to fit an arima model with more than one seasonality (daily, weekly) using a dummy for the weekly seasonality. I found a very good post explaining a similar problem but after many trials I'm running out of ideas with the following error:               Error in optim(init[mask], armaCSS, method = optim.method, hessian = FALSE,  :        non-finite value supplied by optim      My code is the following:               library(lubridate)     start=dmy_hms("25/02/2011 00:00:00")     index=start + c(0:1000) * hours(1)     week = wday(index)     xreg_w = model.matrix(~0+as.factor(week))     colnames(xreg_w) = c("Mon","Tue","Wed","Thu","Fri","Sat","Sun")          freq=24     set.seed(1234)     y=ts(log(35+10*rnorm(1001)),f=freq)          library(forecast)     model = Arima(y, order=c(0,0,0), seasonal=list(order=c(1,0,0), period=24),                    xreg=xreg_w)


title: How to test a large sample if it is normaly distributed?

text: I am dealing with normaly testing of large samples. As stated here: Is normality testing 'essentially useless'? Normally testing is essentially useless, if the sample is too large. Even visual testing cannot give a clear statement about the distribution? So, how to test a large sample if the distribution is normal distributed?


title: Multilevel bayesian with AR1 correlation structure

text: How do I fit a bayesian multilevel model with with AR(1) correlation structure? I am trying to teach myself bayesian modelling and I am wondering how you could specify a multilevel model with an AR(1) correlation structure. e.g how do i get the equivalent from nlme using either `JAGS` or `stan`?               lme(y ~ time + x, random=~time|subject, corr=corAR1())


title: How to calculate mean and standard deviation of a count variable when the raw data is based on frequency categories?

text: I asked 312 people how many times they visited their preferred local supermarket in a month The results are as follows:   * 5% did not visit at all   * 7% visited once a month   * 33% visited twice a month   * 22% visited three times a month   * 15% visited four times a month   * 18% visited five and more times a month In the absence of the actual number of visits (I only have percentage of patrons as above), how do you calculate the mean and standard deviation for reporting purposes.


title: Window models in stream data processing

text: Reading about data stream clustering I met the next terms:   * landmark window model,   * sliding window model,   * damped window. As to sliding window it's clear - oldest data escape the scope, the new data go inside. But what are the consepts of the other two? I can suppose that damped window it's like a buffer, that is flushed after filling but some explanation, that I had found, states that dumped uses decreasing weights of data as function of time.


title: Topic stability in topic models

text: I am working on a project where I want to extract some information about the content of a series of open-ended essays. In this particular project, 148 people wrote essays about a hypothetical student organization as part of a larger experiment. Although in my field (social psychology), the typical way to analyze these data would be to code the essays by hand, I'd like to do this quantitatively, since hand-coding is both labor-intensive and a bit too subjective for my taste. During my investigations about ways to quantitatively analyze free response data, I stumbled upon an approach called topic modelling (or Latent Dirichlet Allocation, or LDA). Topic modeling takes a bag-of-words representation of your data (a term-document matrix) and uses information about the word co- occurrences to extract the latent topics of the data. This approach seems perfect for my application. Unfortunately, when I've applied topic modeling to my data, I've discovered two issues:   1. The topics uncovered by topic modelling are sometimes hard to interpret   2. When I re-run my topic models with a different random seed, the topics seem to change dramatically Issue 2 in particular concerns me. Therefore, I have a two related questions:   1. Is there anything I can do in the LDA procedure to optimize my model fit procedure for interpretability and stability? Personally, I don't care as much about finding the model with the lowest perplexity and / or best model fit -- I mainly want to use this procedure to help me understand and characterize what the participants in this study wrote in their essays. However, I certainly do not want my results to be an artifact of the random seed!   2. Related to the above question, are there any standards for how much data you need to do an LDA? Most of the papers I've seen that have used this method analyze large corpora (e.g., an archive of all Science papers from the past 20 years), but, since I'm using experimental data, my corpus of documents is much smaller. I have posted the essay data here for anyone who wants to get his or her hands dirty, and I have pasted the R code I'm using below.               require(tm)     require(topicmodels)          # Create a corpus from the essay      c <- Corpus(DataframeSource(essays))     inspect(c)          # Remove punctuation and put the words in lower case     c <- tm_map(c, removePunctuation)     c <- tm_map(c, tolower)          # Create a DocumentTermMatrix.  The stopwords are the LIWC function word categories     # I have a copy of the LIWC dictionary, but if you want to do a similar analysis,     # use the default stop words in tm     dtm <- DocumentTermMatrix(c, control = list(stopwords =        c(dict$funct, dict$pronoun, dict$ppron, dict$i, dict$we, dict$you, dict$shehe,              dict$they, dict$inpers, dict$article, dict$aux)))          # Term frequency inverse-document frequency to select the desired words     term_tfidf <- tapply(dtm$v/rowSums(as.matrix(dtm))[dtm$i], dtm$j, mean) * log2(nDocs(dtm)/colSums(as.matrix(dtm)))     summary(term_tfidf)          dtm <- dtm[, term_tfidf >= 0.04]          lda <- LDA(dtm, k = 5, seed = 532)     perplexity(lda)     (terms <- terms(lda, 10))     (topics <- topics(lda))      **Edit:** I tried modifying `nstart` as suggested by Flounderer in the comments. Unfortunately, as shown below, even setting `nstart` to 1000 results in topics that vary quite dramatically from random seed to random seed. Just to emphasize again, the only thing I'm changing in the estimation of the two models below is the random seed used to start model estimation, and yet the topics do not seem to be at all consistent in these two runs.               lda <- LDA(dtm, k = 5, seed = 535, control = list(nstart = 1000))     (terms <- terms(lda, 10))                Topic 1         Topic 2      Topic 3      Topic 4       Topic 5            [1,] "international" "ethnicity"  "free"       "credit"      "kind"             [2,] "communicate"   "true"       "team"       "mandatory"   "bridge"           [3,] "gain"          "asians"     "cooperate"  "music"       "close"            [4,] "use"           "hand"       "order"      "seen"        "deal"             [5,] "big"           "hold"       "play"       "barrier"     "designed"         [6,] "communication" "effective"  "big"        "stereotypes" "effort"           [7,] "america"       "emphasis"   "beginning"  "asians"      "implemented"      [8,] "chinese"       "halls"      "china"      "fantastic"   "websites"         [9,] "ethnicity"     "minorities" "difference" "focusing"    "planned"         [10,] "networks"      "population" "easier"     "force"       "body"          lda <- LDA(dtm, k = 5, seed = 536, control = list(nstart = 1000))     (terms <- terms(lda, 10))                Topic 1       Topic 2         Topic 3        Topic 4       Topic 5          [1,] "kind"        "international" "issue"        "willing"     "play"           [2,] "easier"      "ethnicity"     "close"        "use"         "trying"         [3,] "gain"        "communication" "currently"    "hand"        "unity"          [4,] "websites"    "communicate"   "implemented"  "networks"    "decision"       [5,] "credit"      "bridge"        "particularly" "stereotypes" "gap"            [6,] "effort"      "america"       "credit"       "communicate" "normally"       [7,] "barriers"    "connection"    "fulfill"      "came"        "asians"         [8,] "effects"     "kind"          "grew"         "asians"      "created"        [9,] "established" "order"         "perspectives" "big"         "effective"     [10,] "strangers"   "skills"        "big"          "budget"      "prejudice"


title: Significance of (GAM) regression coefficients when model likelihood is not significantly higher than null

text: I am running a GAM-based regression using the R package gamlss and assuming a zero-inflated beta distribution of the data. I have only a single explanatory variable in my model, so it's basically: `mymodel = gamlss(response ~ input, family=BEZI)`. The algorithm gives me the coefficient $k$ for the impact of the explanatory variable into the mean ($\mu$) and the associated p-value for $k(\text{input})=0$, something like:               Mu link function:  logit                                                    Mu Coefficients:                                                                         Estimate  Std. Error  t value   Pr(>|t|)                       (Intercept)  -2.58051     0.03766  -68.521  0.000e+00                       input        -0.09134     0.01683   -5.428  6.118e-08      As you can see in the above example, the hypothesis of $k(\text{input})=0$ is rejected with high confidence. I then run the null model: `null = gamlss(response ~ 1, family=BEZI)` and compare the likelihoods using a likelihood-ratio test:               p=1-pchisq(-2*(logLik(null)[1]-logLik(mymodel)[1]), df(mymodel)-df(null)).      In a number of cases, I get $p>0.05$ even when the coefficients at input are reported to be highly significant (as above). I find this quite unusual -- at least it never happened in my experience with linear or logistic regression (in fact, this also never happened when I was using zero-adjusted gamma with gamlss). My question is: can I still trust the dependence between response and input when this is the case? Many thanks!


title: QR factorization: floating-point issues

text: I've implemented QR factorization based on Householder reflections (for the purposes of computing the OLS fit). Mathematically, the $R$ matrix is upper triangular. However, due to floating-point issues I typically end up with small non-zero entries below the diagonal. What should I do with them - leave alone or set to zero, and why? For example, the following input:               [[ 1  1 10]      [ 3  3  5]      [ 2  4 -6]      [ 1 10  8]]      produces the following $R$ matrix:               [[ -3.87298335e+00  -7.22956891e+00  -5.42217668e+00]      [ -8.88178420e-16   8.58681159e+00   4.86793026e+00]      [ -4.44089210e-16  -3.33066907e-16   1.31111882e+01]]      Note the non-zeros below the main diagonal. For comparison, both `numpy` and `R` return exact zeros, so I am wondering if there are reasons for me to do the same.


title: Two sample Kolmogorov-Smirnov test and p-value interpretation

text: I have a question regarding the interpretation of resulting p-values of a two sample Kolmogorov Smirnov test. Basis of my analysis is to try to identify groups that show a difference in their distribution difference compared to totality. I used a two sample Kologorov Smirnov Test in R to do so. Sample sizes:               Full = 2409     Group_1 = 25     Group_2 = 26     Group_3 = 33     Group_4 = 43      Dataset plots: ![Boxplots of dataset distributions](http://i.stack.imgur.com/rQqU0.png) "Other" contains a collection of groups containing less than 20 datapoints. The resulting p-values I get when I compare each "Group" against "Full - Group" are the following:               Group 1: 2.6155e-002     Group 2: 2.1126e-003     Group 3: 7.2113e-002     Group 4: 7.3466e-003      How can I interpret these results - especially with regards to the low number of datapoints per group as well as the difference in sample size for Full (N=2409) and Groups (N=25-43)? Is the choice of a KS test good or might another test be more appropriate in this case?


title: How to determine the thresholds in oredered logistic regression

text: from what I have studied in the data mining course (please correct me if I'm wrong) - in logistic regression, when the response variable is binary, then from the ROC curve we can determine the threshold. Now I'm trying to apply the logistic regression for an ordinal categorical response variable with more than two categories (4). I used the function `polr` in r:               > polr1<-polr(Category~Division+ST.Density,data=Versions.data)     > summary(polr1)          Re-fitting to get Hessian          Call:     polr(formula = Category ~ Division + ST.Density, data = Versions.data)          Coefficients:                    Value Std. Error t value     DivisionAP   -0.8237     0.5195  -1.586     DivisionAT   -0.8989     0.5060  -1.776     DivisionBC   -1.5395     0.5712  -2.695     DivisionCA   -1.8102     0.5240  -3.455     DivisionEM   -0.5580     0.4607  -1.211     DivisionNA   -1.7568     0.4704  -3.734     ST.Density    0.3444     0.0750   4.592          Intercepts:         Value   Std. Error t value     1|2 -1.3581  0.4387    -3.0957     2|3 -0.5624  0.4328    -1.2994     3|4  1.2661  0.4390     2.8839          Residual Deviance: 707.8457      AIC: 727.8457        How should I interpret the Intercepts? and how can I determine the threshold for each group? Thanks


title: Excel spreadsheet for ordinal logistic regression

text: Does anyone know of an existing implementation of ordinal logistic regression in Excel?


title: One-class SVM vs. exemplar SVM

text: I understand that one-class SVMs (OSVMs) were proposed with the absence of negative data in mind and that they seek to find decision boundaries that separate a positive set and some negative anchor point, say the origin. A recent work proposes Exemplar SVMs (ESVMs) which trains a "single per- category classifier" which is claims to be different from OSVMs, in that ESVMs do not "require mapping the exemplars into a common feature space over which a similarity kernel can be computed". I do not quite understand what this means and how ESVMs differ from OSVMS. And so how do they differ? And how this is this similarity kernel computation avoided in ESVMs?


title: $\nu$-svm parameter selection

text: For the $\nu$-SVM (for both classification and regression cases) the $\nu \in (0;1)$ should be selected. The LIBSVM guide suggests to use grid search for identifying the optimal value of the $C$ parameter for $C$-SVM, it also recommends to try following values $C = 2^{-5}, 2^{-3}, \dots, 2^{15}$. So the question is, are there any recommendations for values of the $\nu$ parameter in case of $\nu$-SVMs?


title: Anderson-Darling test - what is adjustment good for

text: I am trying to understand and implement Anderson-Darling test for normal distribution in C++. I've found Octave/Matlab implementation of Anderson- Darling test here, but the critical values differ from values listed at wikipedia. Authors of mentioned source code had made some kind of adjustment, before they tested the hypothesis. I would like to know for what is this adjustment good for, why is it made?


title: How to calculate odds ratio and p-values for Fisher's test?

text: I'm trying to perform Fisher's exact test in R. I'm calculating odds ratio like this $$\frac{(\text{no. of successes in my set}) \cdot (\text{no. of failures in background set})}{\text{(no. of successes in background set)} \cdot (\text{no. of failures in my set})}$$ In the denominator, I get 0 for some cases, because there are 0 failures in my set for some cases. Can anyone tell me how to calculate p-values from odds ratio? Also what does it mean if we say odds ratio should be greater or less than 1? If I want to test over-representation in my set against the background set what should be the odds ratio and how can I calculate that?


title: Length weighted median

text: In the field of genome assemblies in biology (the problem being reconstructing the genome from many, short random pieces of it, where a genome is one or more long strings of a limited alphabet), there is a metric called 'N50' for a set of strings (pieces of DNA sequence). The N50 length is defined as 'the length of the string for which 50% of all characters are in strings of that length or longer'. See for example http://seqanswers.com/forums/showthread.php?t=2332. I have sometimes seen this metric referred to as the 'length weighted median'. Is this correct?


title: Trying to test the Significance of the Seattle Seahawks Home Field Advantage

text: Basically I have an interest in sports and as part of that do a Microsoft Excel based version of Kenneth Massey 1997 thesis (http://masseyratings.com/theory/massey97.pdf) method of ranking sports teams, that also allows for a Home Field Advantage (HFA) rating (so I have Solver within Excel solve the least squares method that I have set up based on game results to derive a best fit for the means or ratings for the teams and a common Home Field Advantage). Due to a post on a forum elsewhere, I decided that I wanted to try and put a number on the Seattle Seahawks (NFL) team’s Home Field Advantage (HFA), which allegedly is higher than the league average. Where as my model usually allows for all of the seasons games to solve the least squares problem (to develop an accurate rating), I took data over the past 6 seasons, kept my end of season ratings constant from each of those seasons (derived on results from over 250 games - 16 games for each team) and got Excel to solve the individual Home Field Advantages based on just each of the teams home games (8 per season). The individual HFA's then averaged out to the league average solved by my original model. So effectively I fitted my problem to the model (I think?). I used the Swinburne University of Technology by Paper by Stephen Clarke into “Home Field Advantage in the Australian Football League” (http://researchbank.swinburne.edu.au/vital/access/manager/Repository/swin:1561?exact=creator%3A%22Clarke%2C+Stephen+R.%22) as a basis for the concept and in particular it’s reference to Harville, D. A., and Smith, M. H. (1994). “The home-court advantage: How large is it, and does it vary from team to team?” The American Statistician, 48, 22-28. In short taking data over 6 seasons (2007/2008 – through to 2012/2013), I have found 4 teams that are over two points different to the league average (and 28 that were within 2 points which I didn’t think was a big enough difference) and these teams seem to on average be in either the top 10 or bottom 10 of Home Field Advantages; At the top; Seattle Seahawks – average HFA rank 7th – average difference to league average HFA of 2.58: +3.47 points (e.g. Seattle’s HFA = 6.05) Baltimore Ravens – average HFA rank 9th – average difference to league average HFA: +2.64 points And at the bottom; Carolina Panthers – average HFA 23rd – average difference to league average HFA: -2.11 points Miami Dolphins – average HFA 24th – average difference league average HFA: -2.61points (so on average HFA is marginally negative and the only team that this would affect). My query is; is there a way to test the significance of the 3.47 points above the league average value of the Seattle Seahawks individual HFA? Is it significant? Is the data set of games big enough (6 seasons – 48 games per team)? Am I right to disregard anything below 2 points off of the average as not being significant enough (that it might just be caused by a one season blip)? I am also not looking for the causation of the HFA (which I am expecting to be a combination of familiarity as the stadium has been open only ten years and travel as Seattle is the most isolated of NFL teams). I just want to prove that it is significant (and then I can go after the causes). I don’t have a great maths background (just an interest in sports) so I would be grateful for any advice as I think being consistently 3 points above the league average must be a significant advantage. Thanks,


title: Density Clustering

text: I am looking for a clustering algorithm. My idealized dataset looks like this:   ![http://i.stack.imgur.com/bSlNU.png](http://i.stack.imgur.com/pqEpC.png) The clustering result should look like the Rapidminer density plot:   ![http://i.stack.imgur.com/Sk3xK.png](http://i.stack.imgur.com/uMRQ4.png) Means 3 or 4 clusters should be the clustering result: One Cluster (1|1) to (5|6), one or two for the points on Z = 5 and one for the overlaying area (6|1) to (11|6). Based on "density" I gave the `DBSCAN` R package a shot. The best result I could get was the one in the above picture. So not exactly what I expected, as the overlaying area was not recognized as a separate cluster. Any ideas which algorithm would provide the results I expect? ![enter image description here](http://i.stack.imgur.com/pqEpC.png)


title: What is the difference between GLM and GEE?

text: Whats the difference between a GLM model (logistic regression) with a binary response variable which includes subject and time as covariates and the analogous GEE model which takes into account correlation between measurements at multiple time points? My GLM looks like:               Y(binary) ~ A + B1X1(subject id) + B2X2(time)                    + B3X3(interesting continuous covariate)      with logit link function. I'm looking for a simple (aimed at the social scientist) explanation of how and why time is treated differently in the two models and what the implications would be for interpretation.


title: Interpretation of Conditional Density Plots

text: I would like to know how to correctly interpret conditional density plots. I have inserted two below that I created in R with `cdplot`. For example, is the probability of the _Result_ being equal to 1 when _Var 1_ is 150 approximately 80%? ![conditional density plot](http://i.stack.imgur.com/hst5c.png) The dark grey area is that which is the conditional probability of the `Result` being equal to 1, right? ![conditional density plot](http://i.stack.imgur.com/y3KbF.png) From the `cdplot` documentation: > cdplot computes the conditional densities of x given the levels of y > weighted by the marginal distribution of y. The densities are derived > cumulatively over the levels of y. How does this accumulation affect how these plots are interpreted?


title: Entropy and information content

text: 1. I am curious to know about the relation between Entropy and information content of a signal or trajectory of time series. When a system is at equilibrium, then it has maximum entropy.  Does entropy mean a measure of information loss ? Higher entropy == loss in information? In my understanding, Shannon's information theory entropy measures how uncertain we are about the next incoming symbol. Hence, Entropy = measure of uncertainty. Entropy is maximum when we are absolutely unsure of the outcome of the events. So, would this imply that Max Entropy = minimum information?   1. if we consider a trajectory in phase space and each element of the trajectory to be a symbol, then would the evolution of the trajectory mean that information is being created? How do we infer entropy for a trajectory? Thank you for insights and help in advance.


title: Composite variables in SEM

text: I am doing an analysis of latent constructs. I was wondering whether it was acceptable to use composite variables given that the measures associated with each of the latent constructs uses up to 80 items?


title: Number of trees using cforest in caret

text: I have read some of the previous posts about this issue and I still do not understand how caret decides on the number of trees when you use the `model=cforest` from party package. I know that at certain threshold there is no use of building more trees, but how is this optimal number calculated? Is it a default?


title: Calculating expected profits from probabilities.

text: A drug company is about to start clinical trials of two new drugs that it has developed to help people to quit smoking.   Based on experience and current laboratory results, the drug company believes that there is a 70% chance that both drugs will be accepted and a 20% chance only one will be accepted. If both are accepted the drug company will make \$10 million. If one is accepted the drug company will make a \$4 million profit, while if both are rejected they will lose \$2 million. What is the drug company's expected profit? Is it:- $7.2 million? $7.6 million? $10.0 million? None of the above?