title: Deriving the change in a Poisson lambda over time

text: I have two sample data sets which are each a 1x95 vector containing the $\lambda$ of a Poisson. Each slot contains the $\lambda$ for the minute that the slot represents. So t1[1] is the starting lambda, t1[2] is the $\lambda$ after one minute, etc. The value of the lambda slowly decays over the time series. The series only goes from 1 to 95 minutes. Based on these data sets, is it possible to derive the formula for this decay? I would like to have a formula that could then take a value between say 1.1 and 6.0 and return a vector of $\lambda$'s for each slot 1:95. In other words, I would like to have a function :-               new_lambda <- function(beginning_lambda, current_minute)      The first data set is for initial $\lambda$, 1.4, and the second is for $\lambda$, 2.1. (The first element of each vector is the initial $\lambda$, and each subsequence number is the decayed lambda after slot number minutes)               c(1.4, 1.38731693971954, 1.37461497531924, 1.36189387488779,      1.34915340111498, 1.33639331110529, 1.32361335618285, 1.31081328168731,      1.29799282676014, 1.2851517241207, 1.27228969983157, 1.25940647305241,      1.24650175578168, 1.23357525258547, 1.22062666031263, 1.20765566779532,      1.19466195553401, 1.18164519536608, 1.16860505011665, 1.1555411732307,      1.14245320838518, 1.12934078907955, 1.11620353820353, 1.1030410675802,      1.08985297748287, 1.07663885612371, 1.06339827911211, 1.05013080888051,      1.03683599407519, 1.02351336890948, 1.0101624524762, 0.996782748016363,      0.983373742140467, 0.969934903998559, 0.956465684394797, 0.942965514841851,      0.929433806549955, 0.915869949344926, 0.902273310508843, 0.888643233536399,      0.874979036799168, 0.861280012109178, 0.847545423172193, 0.833774503919968,      0.819966456709522, 0.806120450375986, 0.792235618123945, 0.77831105524032,      0.764345816609622, 0.750338914009939, 0.736289313165104, 0.72219593052513,      0.708057629743121, 0.693873217812284, 0.679641440821358, 0.665360979280522,      0.651030442962453, 0.636648365194487, 0.622213196527442, 0.607723297694309,      0.593176931757131, 0.578572255322537, 0.563907308684732, 0.549180004728389,      0.534388116391651, 0.519529262449729, 0.504600891330413, 0.489600262611446,      0.474524425772686, 0.459370195678448, 0.444134124141028, 0.428812466756316,      0.413401143994426, 0.397895695255455, 0.382291224238561, 0.366582333486939,      0.350763045311143, 0.334826705383382, 0.318765864022086, 0.302572128373676,      0.286235976070479, 0.269746517053042, 0.253091184349661, 0.236255325435227,      0.219221651093633, 0.201969474329622, 0.18447362979703, 0.166702888088314,      0.148617533565322, 0.130165475875085, 0.11127559850736, 0.0918453826079598,      0.0717150018886348, 0.0506024484826778, 0.0278795484294726)          c(2.1, 2.08097540957931, 2.06192246297887, 2.04284081233169,      2.02373010167247, 2.00458996665794, 1.98542003427428, 1.96621992253097,      1.94698924014021, 1.92772758618105, 1.90843454974735, 1.88910970957861,      1.86975263367252, 1.85036287887821, 1.83093999046895, 1.81148350169297,      1.79199293330102, 1.77246779304913, 1.75290757517497, 1.73331175984606,      1.71367981257777, 1.69401118361933, 1.67430530730529, 1.6545616013703,      1.63477946622431, 1.61495828418557, 1.59509741866817, 1.57519621332076,      1.55525399111279, 1.53527005336422, 1.5152436787143, 1.49517412202454,      1.4750606132107, 1.45490235599784, 1.4346985265922, 1.41444827226278,      1.39415070982493, 1.37380492401739, 1.35340996576327, 1.3329648503046,      1.31246855519875, 1.29192001816377, 1.27131813475829, 1.25066175587995,      1.22994968506428, 1.20918067556398, 1.18835342718592, 1.16746658286048,      1.14651872491443, 1.12550837101491, 1.10443396974766, 1.08329389578769,      1.06208644461468, 1.04080982671843, 1.01946216123204, 0.998041468920782,      0.97654566444368, 0.95497254779173, 0.933319794791163, 0.911584946541464,      0.889765397635697, 0.867858382983806, 0.845860963027097, 0.823770007092583,      0.801582174587476, 0.779293893674594, 0.75690133699562, 0.734400393917169,      0.711786638659029, 0.689055293517672, 0.666201186211542, 0.643218700134474,      0.620101715991639, 0.596843542883182, 0.573436836357841, 0.549873500230409,      0.526144567966714, 0.502240058075073, 0.478148796033129, 0.453858192560515,      0.429353964105719, 0.404619775579563, 0.379636776524491, 0.35438298815284,      0.32883247664045, 0.302954211494433, 0.276710444695545, 0.250054332132471,      0.222926300347984, 0.195248213812627, 0.16691339776104, 0.13776807391194,      0.107572502832952, 0.0759036727240167, 0.041819322644209)      Any pointers or hints would be greatly appreciated. Thanks.


title: Is Benjamini-Hochberg correction more conservative as the number of comparisons increases?

text: How conservative is Benjamini-Hochberg multiple testing correction relative to the total number of comparisons? For example, if I have a list of 18,000 features for two groups and I perform a Wilcoxon test to get a p-value. I adjust that p-value using Benjamini-Hochberg and next to nothing comes out as significant. I know that Bonferroni correction can be quite conservative as the number of comparisons increases, does Benjamini-Hochberg have the same property?


title: What is the outline for the procedure of model selection, with different models based on likelihood functions?

text: I have two basic models $M_1$ and $M_2$. They each have a likelihood function; $L_{M_1} = f(\mathbf{X}|\mathbf{\theta_1})$ and $L_{M_2} = f(\mathbf{X}|\mathbf{\theta_2})$ (here $\mathbf{X}$ is the data set, and $\mathbf{\theta_i}$ are the parameters for the model). The models are nested since model $M_1$ parameter's are all present in model $M_2$. Specifically, $\mathbf{\theta_1}$ has 2 parameters and $\mathbf{\theta_2}$ has 3 parameters. QUESTION 1)The model $M_1$ is termed a 'null model' in this case as it is the more simple in the nested case? QUESTION 2)I want to do model selection. The more complicated model, $M_2$, has a greater likelihood but I would like it to be statistically justified. Is it correct to assume that the **likelihood-ratio-test** is the suitable choice to apply? QUESTION 3)For the likelihood ratio test, I am supposed to use the parameter values which maximize the likelihood. In the case of model 1 $M_1$ I can do this analytically, but for model 2 $M_2$ I must approximate it via sampling. Is it suitable that I use $\mathbf{\theta_1}^{max}$ and $\mathbf{\theta_1}^{approx.max}$ in the same mechanism for comparison? I applied the likelihood-ratio-test, and get this result (in logs):               LogLikModel1 = -61159.4324     LogLikModel2 = -30577.2917     LogLikRationTest_D = 61164.2813      QUESTION 4)To estimate the significance the chi-squared, $\chi^2$ test must be used? QUESTION 5)If the answer to the question 4) is yes, then must I use this table (link to chisquared and p-values) for mapping the degrees of freedom and chi values to p-values? Are the numbers that the likelihood-ratio-test produce chi values? How do I use these values in logs with the table for chi-squared p-values because the probabilities are too low to not use logs?


title: How to determine how many unique visitors should be in the control group

text: Let's say I want to perform an AB-Test with the following properties:   * We assign unique visitors to each group, not the requests (so each trial corresponds to a unique visitor as long as e.g. the cookie is not deleted)   * We want to measure "items viewed" per unique visitor (side note) The question now is: How to define the percentage of all trials / unique visitors which should be assigned to the control group A so that ...   1. _A_ should be as small as possible   2. _A_ should not be tested significantly different to _B_ according to a suitable test assuming that there is no difference between _A_ and _B_


title: Learning distributions from "tree based" matrix arrangements

text: Disclaimer : I am not a math/stats major. Please pardon & feel free to correct non-standard words I may use. Following is an arrangement I currently have with me, I have a matrix r x c where value of each cell is a binary 0/1. This r x c matrix is "connected" to further r matrices. Such a connection is done by connecting each row of the r x c matrix to another matrix. Each of these matrices is similarly connected to other matrices. One can imagine of a graph where each node is a matrix and each edge is a connection between such matrices. Question : Is there a way to exploit the structure of such a graph to form a relation between the values of the matrices (preferably a graph itself!) ?


title: How to detect changes in amplitude?

text: I have timeseries like this: ![enter image description here](http://i.stack.imgur.com/2oNmT.png) as you can see there are changes regarding the amplitude. Is there a test to check this kind of changes? Important annotations:   1. I do not know if the series have changes in amplitude   2. If there is a change in amplitude I do not know the point of the change   3. The changes can be more then ONE (but I only need to know if there is a change in amplitude, for my tests the numbers of changes is not important)   4. As you can see the means are common   5. I do not have groups I have series in a numeric vector (R vector), I only subdivided the above series in three groups to show the three changes in amplidute, but it is very obvious.


title: Sidak or Bonferroni?

text: I am using a generalized linear model in SPSS to look at the differences in average number of caterpillars (non-normal, using Tweedie distribution) on 16 different species of plants. I want to run multiple comparisons but I'm not sure if I should use a Sidak or Bonferroni correction test. What is the difference between the two tests? Is one better than the other?


title: Poisson or quasi poisson in a regression with count data and overdispersion?

text: I have count data (demand/offer analysis with counting number of customers, depending on - possibly - many factors). I tried a linear regression with normal errors, but my QQ-plot is not really good. I tried a log transformation of the answer: once again, bad QQ-plot. So now, I'm trying a regression with Poisson Errors. With a model with all significant variables, I get:               Null deviance: 12593.2  on 53  degrees of freedom     Residual deviance:  1161.3  on 37  degrees of freedom     AIC: 1573.7          Number of Fisher Scoring iterations: 5      Residual deviance is larger than residual degrees of freedom: I have overdispersion. How can I know if I need to use quasipoisson? What's the goal of quasipoisson in this case? I read this advise in "The R Book" by Crawley, but I don't see the point nor a large improvement in my case.


title: Notation for Random Bernoulli-Like Vector With Fixed Sum

text: I draw a random vector of dimensionality $k$, each dichotomous element of which taking on a value in $\\{0,1\\}$. The probability that any element will be $1$ is captured in the $k$-dimensional parameter vector $\bf{p}$, the elements of which are in $[0,1]$ and sum to $1$ (or are normalized as such prior to use). A special condition is that the random vector must sum to some value $c$. If you're familiar with R, I can obtain such random vectors in R for some $k$, $c$, and $\bf{p}$ by               v <- vector(length=k)     v[sample(1:k, c, prob=p)] <- 1     v      I am hoping someone can point me to a distribution that captures this. Ultimately, I need a **notation** to use to describe this succinctly in a paper. Thank you very much for your time.


title: With regards to asymettric error bars, what is the correct alternative to S.E.M. when data are strictly positive and near zero?

text: I have some data that measures how a substance decays over time. At each time point I have 4 measurements. At time points where there is a lot of substance I use the standard error for error bars and this seems standard practice. However, when the amount of substance nears zero and I have a lot of variability in the measurement the standard error indicates that the lower error bars are values less than zero and this is physically impossible for this substance. When this is the case, how should I calculate where my error bars go? I suppose I'm assuming the underlying distribution of the data is gaussian when I use standard error, but clearly the data is not gaussian when it is near zero, because I can't have negative amounts of substance. Is there another distribution I should be using? Thanks.


title: What is the closest approximation to computing the $X$th percentile (say 95th) from bins of data?

text: I have a few bins of data. From each bin(which has a few thousand data points), I can compute the total number of data points, the mean of the data points, etc.. I can probably compute any/most statistical entities(mean, median, etc.) of each bin(computationally easy). However, I want to find the 95th percentile of all the data points. The total number of buckets are not that many(~90), so again, I have some room for computation amongst individual statistical entities of each bin. For finding the 95th percentile, I would generally sort all the data points and pick out the one at the 95th percentile. That won't work in my case(too many data points, computationally intensive). Is there any close approximation to the 95th percentile that I can compute? There are two phases in my problem - one when I can compute statistical entities for each bin, and the second when I do something with the statistical entities of each bin. Phase 1 and phase 2 happen in different locations, and I can't pass too much data across the two of them. This is turning out to be a mix of a statistical and computer systems design problem. Please let me know if there is a better place for such a problem. :/


title: How do I account for both intra- and inter-individual variance when comparing the control with a conditioned test stimulus?

text: I have data that consists of 10 repetitions of a measurement per person, for all repetitions there is one control value and one test value (after conditioning). This data will be reproduced in multiple people. I am now thinking how to introduce the intra-individual variance into the statistical test. Without repetitions I would simply use a paired t-test if I can expect a normal distribution of the data. I would be very thankful if someone could point me to some literature or explain to me how to analyze this sort of experiment. ### Edit: Answers to Michelle's questions This is a neurophysiological study. The goal is to observe the influence of conditioning peripheral electrical stimulation on the amplitude of reflex responses. Thus there is always a pair of reflexes first with 10 s no activity/stimulation before (control) and afterwards again 10 s no activity/stimulation followed by 10 s of conditioning. Here is a small schematic (. no stimulation, I reflex stimulation and - conditioning stimulation): .........I..........----------I This paradigm is repeated per person 10 times, i.e. 20 measures per person, always one control matched to one test measure. The time lag between the end of conditioning and test is important and will be investigated in the future but for this experiment there is only a fixed delay. Nonetheless I would be thankful to know how to handle this too. I do not believe that there is a known mechanistic process which affects the values.


title: Using collaborative filtering to "clean" data and the other way around

text: I am considering two types of systems - which might have more appropriate names: **Recommender systems:** These recommender systems are based on collaborative filtering methods, both model-based and memory-based. They take a partially filled user/item matrix of ratings of items by users, and try to predict the missing ratings. **Trust systems:** These systems take a filled user/item matrix of ratings, and try to detect a trustworthiness for each user - which can then be used to give more fair ratings. These systems could for example be used to eliminate overly negative users or fake ratings. There seems to be a lot of overlap between these two types of systems. My question is, if it is possible to simulate one system using the other. For example - to simulate a trust system with a recommender system - one could remove a rating $r$, let the recommender system predict it, and see how big a difference the actual and predicted value have. Is it possible to use such a trick? What is the correct terminology for these systems? I am most interested in systems of the second type, but have so far stumbled upon most research papers on the first type.


title: How to specify a lognormal distribution in the glm family argument in R?

text: Simple question: How to specify a lognormal distribution in the GLM family argument in R? I could not find how this can be achieved. Why is lognormal (or exponential) not an option in the family argument? Somewhere in the R-Archives I read that one simply has to use the log-link for the family set to gaussian in the GLM, in order to specify a lognormal. However, this is nonsense because this will fit a non-linear regression and R starts asking for starting values. Is anybody aware how to set a lognormal (or exponential) distribution for a GLM?


title: Difference between the Stopping Criteria

text: I would like to know the difference between below mentioned stopping criteria used in various gradient descent algorithm   1. $\frac{Prev\\_fun\\_value - curr\\_fun\\_value}{Pre\\_fun\\_value} \le tol$   2. $Prev\\_fun\\_value - curr\\_fun\\_value \le tol*max(1,Prev\\_fun\\_value)$ where $Prev\\_fun\\_value$ is the previous function value before the update and $curr\\_fun\\_value$ is the current function value after the update of optimization variable and $tol$ is the tolerance.


title: Using Poisson instead of transforming the data

text: I note with interest an article which has been suggested to me over here: http://blog.stata.com/tag/poisson-regression/ I have the problems outlined in this blog - namely zeros (in the independent variables). The article suggests the use of Poisson even though there is no count data in the dependent variable. What should I be on the look out for in terms of checking that Poisson is suitable for my data and that the results I am obtaining are meaningful? In previous posts I have already indicated the nature of my data. The dependent variable is investment in millions of pounds, whilst the independent variables include age of firm, and a number of indicators. I have also other regressions to run which include a number of counts in the independent variable which also could be zero. In my case, investments tend to be skewed towards lower amounts, and so is the age of firm. There is an issue with using OLS - the normality assumption of residuals does not hold, given the bias towards smaller investments and younger firms.. also there is evidence of heteroskedasticity (ok I can regress with robust standard errors to account for this). I have tried using mboxcox and the ladder feature to try and identify possible transformations but all efforts till now proved useless. If I use iqr - I note that i have < 2% mild outliers, and no severe outliers. Trust you'd be able to help further on this...


title: Can I compare relative frequency between groups and test significance?

text: This is a biological experiments question. I compare the percentage (%; relative frequency) of mutant cells showing phenotype against wild type cell (wt) with repeated times, for example, three independent trials. I: wt (0.02 = 2 out of 100) vs mutant (0.50 = 50 out 100) II: wt (0.027 = 4 out 150) mutant (0.44 = 40 out 90) III: wt (0.01 = 1 out 90) mutant (0.29 = 35 out 120) Finally I want to represent the frequency (%) of mutant phenotype with some error bars, which is typical biological graphic data representation. I used to test the significance of difference between samples with T-test using % number. (eg. 0.02,0.027, 0.01 vs 0.5,0.44,0.29) However, sometimes sample number size from each single experiment can be slightly different. I am not sure my method is correct. Can anybody give advice for analyzing of this kind of data? Thank you in advance.


title: Metrics for comparing estimated lists to a 'true' list

text: I'm wondering what the best ways to compare (possibly ranked) lists when we know what the true ranking is and also the variable that decides the ranking. Say this is the top 10 of a certain list, we know the ranking, and also the variable from which the ranking is obtained (Votes). **True Ranking:**                                    Player Votes     1           Ablett, Gary GC    28     528        Selwood, Joel GE    27     588           Swan, Dane CW    26     301       Johnson, Steve GE    25     120 Dangerfield, Patrick AD    22     236       Hannebery, Dan SY    21     464    Pendlebury, Scott CW    21     502        Rockliff, Tom BL    21     102       Cotchin, Trent RI    19     285         Jack, Kieren SY    19      Now, say I have two lists, produced by models A and B. The models have been trained on a independent set of data to predict the number of votes obtained by each player on a new dataset (from which the True Ranking is associated with). **Model A output:**               1        Ablett, Gary GC    41     528     Selwood, Joel GE    30     588        Swan, Dane CW    29     211     Griffen, Ryan WB    28     464 Pendlebury, Scott CW    24     502     Rockliff, Tom BL    24     641      Watson, Jobe ES    23     301    Johnson, Steve GE    22     102    Cotchin, Trent RI    21     180      Fyfe, Nathan FR    21      **Model B output:**               1           Ablett, Gary GC 29.34127     588           Swan, Dane CW 25.49142     211        Griffen, Ryan WB 22.50983     464    Pendlebury, Scott CW 19.84517     528        Selwood, Joel GE 18.32023     120 Dangerfield, Patrick AD 16.94963     301       Johnson, Steve GE 16.05056     641         Watson, Jobe ES 15.73885     416      Montagna, Leigh SK 15.35478     339      Liberatore, Tom WB 14.50770      What are the best metrics or loss functions for determining whether model A or B's output is closer to the truth? I'm not sure if it is better to compare the rankings, or to compare the differences in votes for each player between lists. Are there optimal ways to do either? or does it depend on how one chooses to weight things? In this setting, the position of the individuals in the list is probably more important than the number of votes obtained, but I imagine this information can still be used in some way. One of my concerns is that the number of ties in votes between individuals increases dramatically as you go down the 'True list'. In which case I imagine simply taking comparisons on the first 20 or 50 ranked entries might be helpful. Furthermore, the further down the list, the less important it is to be ranked correctly. e.g., It matters a great deal if the items in position 1 and 2 are swapped, but it's basically irrelevant if items 15 and 16 are swapped. I am familiar with rank correlation methods, but are there other measures that might be more appropriate here? Are the measures mentioned here useful for solving this problem? Cheers.


title: JAGS meta-analysis

text: I'm using JAGS to run a meta-analysis and I've run into an issue. I have calculated log-response ratios and errors for about 177 studies, where studies fall into one of 8 groups. I'm interested in estimating an overall effect as well as a effects for each group. I've discovered that if I use the model               model{       for(i in 1:N){         RR[i] ~ dnorm(yhat[i], prec_y[i])         yhat[i] <- mu + gamma[type[i]]       }            for(j in 1:J){         gamma[j] ~ dnorm(gamma_mu, gamma_prec)         gamma_true[j] <- gamma[j] - mean(gamma[])       }            mu_true <- mu + mean(gamma[])            mu ~ dnorm(0, 0.001)       mu_prec ~ dgamma(0.01, 0.01)       gamma_mu ~ dnorm(0, 0.001)       gamma_prec ~ dgamma(0.01, 0.01)     }      than the model doesn't fit my data well, because the estimate for Plant-Mych is just way off: ![sample image](http://www.natelemoine.com/not_working.png) Alternatively, if I use the model:               model{       for(i in 1:N){         RR[i] ~ dnorm(yhat[i], prec_y[i])         yhat[i] ~ dnorm(gamma[type[i]], gamma_prec[type[i]])       }            for(j in 1:J){         gamma[j] ~ dnorm(overall_mu, mu_prec)         gamma_prec[j] <- pow(gamma_sigma[j], -2)         gamma_sigma[j] ~ dunif(0, 100)       }            overall_mu ~ dnorm(0, 0.001)       mu_prec ~ dgamma(0.01, 0.01)          }      then it fits much better: ![sample image](http://www.natelemoine.com/working.png) Is the fundamental difference that, in the first bad model, the groups are fixed effects whereas in the second model the groups are random? If so, why does this make such a big difference? Is there a way to code the first model so that it fits better? I appreciate the help!


title: What does statistical power mean when we are interested in the probability of correctly not rejecting the null hypothesis?

text: It seems that statisticians are usually only interested in statistical power. In other words, they are interested in the probability of correctly rejecting the null hypothesis.   * What if we are interested in the probability of correctly not rejecting the null hypothesis?   * Would we change the hypothesis test so that we are computing a power?


title: What is the difference between Inference and Machine Learning?

text: I have seen some classes in my University labeled as an "Inference class" and others as "Machine learning" classes, but I was not sure if appreciated the core difference between these two labelings? I normally think of machine learning as the study of systems/algorithms that are able tp extract information from data. It uses techiniques from a broad list subjects, like statistics, computer science, optimization, etc. However, I was not sure how statistical inference differed from Machine Learning. Thanks for the help in advance! :)


title: Mixed model specification with nlme in R

text: i'm having a hard time finding the right model for my data. The study is about the effect of a drug (more precise: the effect of the dose of a drug) and the time on the subject (rat). So i got the following variables:   * Dependent Variable: Relative Change   * Independent Variables: Dose (Levels: 1, 2, 4, 10) and Time (1,2,3)   * Further we got the subject ID (1 to 6) resulting in data like this:               ID Dose Time Change     1  1    1    0.342     1  1    2    0.212     1  1    3    0.121     1  2    1    0.423     1  2    2    0.321     1  2    3    ...     1  4    1    ...     ....     6  10   1    0.121     6  10   2    0.221     6  10   3    0.111      As you can see the different doses were applied to all subjects. Three measurements were made for each subject+dose after 1, 2 and 3 minutes. My idea for a model: $$ y_{ijt} = \beta_0 + \beta_1 dose_j + \beta_2 time_t + \beta_3 dose_j*time_t + u_{i} + \epsilon_{ijt}$$ with Dose and Time being fixed effects and $u_i$ being the subject-specific random intercept. Of course this model does not respect the repeated measurements. So I would expect the residuals to be correlated (over the three time-points), modeling this with an adequate residual variance structure like the following: $$ \epsilon_{i} \sim N(0, R_i)$$ with $\epsilon_i$ being a vector of all residuals of subject i and $R_i$ being a block-diagonal matrix with diagonal elements $$ R_{i,diag} = \begin{pmatrix} \sigma^2 & \sigma_{1,2} & \sigma_{1,3} \\\ \sigma_{1,2} & \sigma^2 & \sigma_{2,3} \\\ \sigma_{1,3} & \sigma_{2,3}& \sigma^2 \end{pmatrix}$$. $$ R_i = \begin{pmatrix} R_{i, diag} & 0 & 0 & 0 \\\ 0 & R_{i, diag} & 0 & 0 \\\ 0 & 0 & R_{i, diag} &0 \\\ 0 & 0 & 0 & R_{i, diag} \\\ \end{pmatrix}$$ I would like to know:   1. Is this an adequate way to model this situation or is there a more simple way to do that (would be great! :-)).   2. If this is a right way: Is there a way to estimate this using the lme-function in R (nlme-package)? The problem I have: how might I specify this residual variance-covariance-structure? I got the impression that I need to implement my own varFunc, which doesn't seem to be fun. Concluding remark: At the moment I have no access to "Mixed-Effects Models in S (Pinheiro and Bates)". It will arrive next week. :-)


title: PCA vector interpretation

text: I have various object vectors (PC1, PC2, PC3) representing objects ater normalisation and PCA. I also have 'axis' vectors (PC1, PC2, PC3) representing the axis that the objects were originally placed along e.g. hot/cold, like/dislike. What I would like to know is if I do a projection of the object vectors on to the 'axis' vectors does this effectively give a representation of the original (normalised for all datasets) data? Can I then take object- axis projections for many entries to give average and s.d. values? I add that the PCA and normalisation was done with multiple entries for the same bunch of objects and axis i.e. 20 people x 16 objects x 5 description axis (but description axis considered seperately i.e. vectors all different). Unfortunately it was also done in software which is rather black box in what it does (please dont just tell me to use something else though).


title: Variable selection with restricted cubic splines

text: Is there any function in R for doing variable selection (backward elimination) in a multiple logistic regression using restricted cubic splines like **mvrs** procedure for STATA?


title: Heavy tailed multi-variate distribution

text: I have been asked to look at some spatial data, which I have assumed is multi- normally distributed in 3 dimensions (and to a reasonable approximation this appears to be true). Actually, all I've been asked to do is provide a table so that someone can specify a probability and I tell them how many standard deviations of their data they need to have that coverage so I wasn't expecting this to take long (i.e. 1.96sd = 95% in 1D). On a closer inspection of the data it is rather heavy-tailed and based on what I think it's wanted for these are the interesting data so I think I possibly shouldn't be using the normal distribtuion, making this a more interesting task. The only two-tailed, tail heavy distribution I know of is the t-distribution and a quick google search seems to imply that there is not a natural multi-variate extension of this although I could just use the one on wiki. Is there a multi-variate t-distribution that's considered 'best' or is there a better distribution that I should be considering? Thanks.


title: Modelling longitudinal data

text: We have longitudinal data on children(n<20) in which we measure different quantities A,B,C,D (like distance walked, time spent in school etc.). These are all continuous variables. We measure these variables daily and we have been measuring them for more than one and half years. Now, based on these variables (A,B,C,D), I am trying to model individual growth trajectories of each child. Ideally, I want to have a single curve that explains how these variables have evolved with respect to time. For now, I am not concerned if different child have different trajectories. I have been studying multi-level models to model growth curves. These techniques were not helpful particularly because I am not interested in modelling individual variables separately and my time series data is quite long(n>500 for each individual). I want to find a relation between all these variables and time. I have been wandering for quite a time. I am looking for getting started papers or some ideas.


title: Interpretation of very small GEE Interactions

text: I have multiple GEE models predicting learning scores from the interactions of psychophysiological data of subjects, e.g. features of Electroencephalograph (EEG), respiration, etc., in a 2 condition design. Thus for example one model is specified (in SPSS syntax) as follows (covariates are EEG delta power, delta baseline, Respiration rate, and baseline for Resp).               * Generalized Estimating Equations.      * Interaction:norm2_RESP*norm_DELTA.      GENLIN L_Score_Exp BY COND (ORDER=ASCENDING) WITH norBL_RESP norBL_DELTA norm2_RESP norm_DELTA      /MODEL COND norBL_RESP norBL_DELTA norm2_RESP norm_DELTA norm2_RESP*norm_DELTA INTERCEPT=YES DISTRIBUTION=NORMAL LINK=IDENTITY      /CRITERIA SCALE=MLE PCONVERGE=1E-006(ABSOLUTE) SINGULAR=1E-012 ANALYSISTYPE=3(WALD) CILEVEL=95 LIKELIHOOD=FULL      /REPEATED SUBJECT=SUBJ WITHINSUBJECT=TRIAL*MIN SORT=YES CORRTYPE=AR(1) ADJUSTCORR=YES COVB=ROBUST MAXITERATIONS=100 PCONVERGE=1e-006(ABSOLUTE) UPDATECORR=1      /MISSING CLASSMISSING=EXCLUDE      /PRINT CPS DESCRIPTIVES MODELINFO FIT SUMMARY SOLUTION (EXPONENTIATED).      The output of this model is an interaction effect which is significant, has B=-.003, and effects for the interaction components, DELTA and RESP, which are significant with B=.003 and B=.008 respectively. Intercept is 10.557 The Bs are very small relative to the intercept, and the B for the interaction has the opposite sign to each of the Bs for components of the interaction alone. The question is, based on the model values how does one interpret the interaction*? Any insights you can lend would be greatly appreciated! :) *I have interpreted graphically by binning but would like to understand the model itself better.


title: Interval censored Cox proportional hazards model in R

text: Given interval censored survival times, how do I perform an interval censored Cox PH model in `R`? An rseek search turns up the package `intcox`, which no longer exists in the `R` repository. I'm almost positive the `coxph` function in the `survival` package cannot handle interval censored survival data. Also, I don't want to impute the data and then use the `coxph` function. This method underestimates the standard errors of the coefficients because you are ignoring the uncertainty of the interval censoring.


title: Is "discriminant function" a synonym for "classification function"

text: "discriminant function" and "classification function" are two terms used in literature to denote a a function that maps a feature vector into a discrete class variable. I presume "discriminant function" has it's origin in statistics and "classification function" in machine learning. Is there any reason these terms cant be used as synonyms?


title: How do important and insignificant variables impact model?

text: I am working to build a model using logistic regression. There is one variable which has strong predictive power. But based on business rule, I cannot include this variable in the model. On the other hand, there are some variables which are insignificant in the model. But from business decision, they are added in the model. So my question is how the model will be impacted if I remove an important variable? And how will it be impacted if I add some insignificant variables? Is there a way to measure the impact?


title: Is SVM-RFE a "filter" or a "wrapper" feature selection algorithm?

text: I am wondering whether SVM based recursive feature elimination (Guyon et al., 2000) would be considered a filter or a wrapper method. On the one hand, it "wraps" around an SVM to get the lowest-weight features. On the other hand, the trained model is never applied to any test set. What do you think?


title: Aligned rank transform on time series data

text: I have some time-series data that I wish to test to see if a certain factor has an effect. There are few samples per time step (10 on average), but the whole time-series has 306 samples. When divided into groups based upon a factor the result is non-normal and non homogeneous. I wish to test the effect of the factor, and if possible, determine a linear trend over the time-series accounting for this effect. To test the effect, I performed an aligned rank transform on the data, then used an ANOVA. However this does not take into account time. To do this I tried a General Linear Model on the transformed data using the factor as a fixed effect and time as a covariate. Does this sound correct? What is the correct way to include time information in a factorial analysis when the number of samples per unit time is quite low? Edit: Additional information: **Time Series** ![Time Series](http://i.stack.imgur.com/UD51S.png) **Factor** ![Factor](http://i.stack.imgur.com/ILGk4.png) Note: software used is SPSS, stats experience is beginner. **Download Data Here:** Data


title: Meta-analysis in R of pre-post treatment design studies

text: I would like to conduct a meta-analysis in R of studies using a pre-post treatment design. In every individual study 4 different subgroups of patients have undergone a pharmacological treatment and their symptoms are measured on a continuous scale. For every study mean and sd of the effect size (cohens d) for the treatment effect is available for every group (this makes the analysis for complicated then usual meta-analysis in which only a single outcome measure per study is used). What is an appropriate way to investigate the effect of treatment (difference between pre- and post-measurement) as well as the effect of group as well as the group-treatment interaction? Would the below example be a easy solution? Here is some R code for dummy data:               set.seed(123)     a <- c(rep("study1",4), rep("study2",4), rep("study3",4), rep("study4",4),       rep("study5",4), rep("study6",4))     b <- rep(c("group1", "group2", "group3", "group4"),6)     data <- as.data.frame(cbind(a,b))     names(data) <- c("study", "group")     data$d <- rep(c(1,2,3,4),6) + rnorm(24,0,0.5)         data$sd <- rep(1, 24) + rnorm(24,0,0.25)               library(metafor)          model_1 <- rma(d, sd, data=data)     summary(model_1)     forest(model_1)          model_2 <- rma(d, sd, mods=group, data=data)     summary(model_2)     forest(model_2)


title: Experiment design: Can unbalanced dataset be better than balanced?

text: I'm on the stage of experiment design of some biomedical time-course study. Let's say we will have 2 groups of subjects - case and control. The total number of subjects is limited (for example, 30), and we need to decide how it's better to split subjects between case and control. One way is to have balanced dataset 50% of case vs. 50% of control. Another way is to assign 2/3 of subjects to case and 1/3 to control, since the case subjects would be more interesting. What experiemnt design would you suggest? Which issues/pitfalls should be considered in both case? Should the experiment design depend on questions asked?


title: Is residuals autocorrelation always a problem?

text: I read that OLS underestimates variance when residuals are autocorrelated. I see why autocorrelation would be a problem in time series analysis, in the sense that the coefficient are not efficient because we're not including all the potential predictors. But is there a mathematical problem as well? For example, we want to predict the used-cars sales margins. The data set includes each vehicle make, model, mileage p/gallon, price, options etc. and the final sale price. For some reason the catalog has been sorted by car make and year/model, so adjacent observations will likely have similar sale numbers. Is autocorrelation a problem in this case?


title: Testing that two models from two different data sets are independent

text: I have two different sets $(n_1=47$, $n_2=23)$ obtained under different conditions. I fit two different functions $(\text{Fu}_1$, $\text{Fu}_2)$ in MATLAB. $\text{Fu}_1$ was fit using the first data set, and $\text{Fu}_2$ was fit using the second data set. $\text{Fu}_1(x,y) = A_1 + B_1x + B_2y$   and   $\text{Fu}_2(x,y) = A_2 + B_3x + B_4y$ When I plot these functions, $\text{Fu}_1$ with the first data set $(n_1=47)$ and $\text{Fu}_2$ with second data set $(n_2=23)$, they are separated in space. How do I test in MATLAB that the models are indeed different?


title: Log transforming dependent variables

text: Does anyone have a clear explanation or reference as to why it is only important to transform the dependent variable and not necessarily the independent variables for improved skewness and kurtosis values?


title: Twitter user data - Fixed Effects Logit the right/only choice?

text: I'm wondering whether a fixed-effects (logit) model fits the purpose of my study, and would like to get some feedback or ideas which other models could potentially be useful: I have an unbalanced panel (daily observations, 3 months time period) with Twitter user data. Among the independent variables are   * number of followers   * is the user verified (binary)   * number of tweets on that day   * time between tweet   * ... the dependent variable is binary (0,1). As far as I understand, the fixed effects model does **not** consider the variance between the individuals. However, isn't the variance **between** the user what really matters here? Given that I suspect that there is much less variation for the independent variables **within** each user (e.g. the number of followers does not increase so much during the study period. Also whether or not a user is a verified user does not change). Still, I understand the advantages of the fixed-effects model, given that there are likely individual time-constant differences between the users (e.g. intelligence, talent,...) which the FE model would control for. Is there a model that accounts for individual fixed effects but still considers between variance? Also, in case I wanted to do a follow-up cluster analysis with the variables that turned out to be statistically significant in the regression, wouldn't I need to consider the between variance in the regression in the first step, because any cluster algorithm essentially tries to separate ( **between** ) users. Thanks a lot for your help and ideas!


title: Non-fixed seasonality in time series

text: I believe some time series data I am looking at shows seasonality according to general elections (UK, so this can be 4 or 5 years). How do I remove this effect in R?


title: What is SVM regression? Is it for regression or classification?

text: I'm trying to understand what is SVM regression. It's used for classification or regression? Can someone give an intuitive understanding of it?


title: Out-of-sample embedding into principal coordinate space

text: I'm trying to project a point into an existing PCOA (Principal Coordinates Analysis) space (in `R`). I am under the impression this must be possible, but I can't figure out how to go about it. Here's how far I've gotten (a toy example):               x <- c(1:10)     y <- c(10:1)     z <- c(rnorm(10,mean=0,sd=2),rnorm(10,mean=10,sd=2))     m <- cbind(x,y,z)          d <- dist(m)     r <- pcoa(d)     biplot(r,m)      The biplot generates the representation I want. Now, given a new point P=(x,y,z) I would like to project it into the above space. The reason I need this and can't simply add this point to the original matrix is that this new point is going to be an outlier and would change the projection of the original space. What I want is to know where this point ends up relative to the untainted representation. Also note that I don't actually use a Euclidean distance in reality, so doing a PCA is not an option.


title: Metropolis-Hastings within Gibbs sampling

text: Suppose we have the following classical normal linear regression model: $$y_i = \beta_1 x_{1i} + \beta_2x_{2i} + \beta_3x_{3i} + e_i$$ where $e_{i} \sim iid.N(0, \sigma^2)$ for all $i = 1, 2, \cdots, n$ and $x_{1i} = 1$ for all $i = 1, 2, \cdots, n$. I artificially generate data values for both $x_{2i}$ and $x_{3i}$ such that $x_{2i}$ and $x_{3i}$ $\sim$ $iid N(0,1)$ for all $i = 1, 2, \cdots, n$. Define $\boldsymbol{\beta} = (\beta_1, \beta_2, \beta_3)'$ and assume a non- informative prior of the form $p(\boldsymbol{\beta}, \sigma) \propto \frac{1}{\sigma}$. * * * Part 1: (Already solved here but sets the context for the next part) Use the conditional posterior pdfs $p(\boldsymbol{\beta}|\sigma, \mathbf{y})$ and $p(\sigma|\boldsymbol{\beta},\mathbf{y})$ to run a Gibbs sampler and estimate the posterior pdf for $\displaystyle \psi = \frac{\beta_2 + \beta_3}{\sigma^2}$ with $n=20$ as the sample size and $M=1000$ as the number of draws after a burn in period of $B = 100$. Part 2: (Attempted but unsure regarding my solution). Insert a Metropolis Hastings sub-step into the program in part 1 at iteration $j-1$ of the Gibbs algorithm to draw from the conditional posterior of $\boldsymbol{\beta}$. Use the following candidate pdf: $$q(\boldsymbol{\beta}|\sigma^{j-1}, \mathbf{y})=(2\pi)^{-\frac{3}{2}}\left|\sigma^{(j-1)2}\boldsymbol{\Omega}^{-1}\right|^{-\frac{1}{2}}\exp\left[-\frac{1}{2\sigma^{(j-1)2}}(\boldsymbol{\beta}-\boldsymbol{b})'\boldsymbol{\Omega}(\boldsymbol{\beta}-\boldsymbol{b})\right]$$ where $\sigma^{(j-1)}$ is the draw of $\sigma$ at iteration $j-1$ of the Gibbs algorithm. In specifying $\boldsymbol{\Omega}^{-1}$, mimic the basic magnitudes of the values in $(\boldsymbol{X}'\boldsymbol{X})^{-1}$ but making the numbers different. * * * The following is my attempt at part 2 and coded in R. Since the MH is applied within an outer Gibbs algorithm, the MH chain does not need to be iterated as convergence occurs immediately and only one value needs to be draw from the MH algorithm. The steps of the MH (within the Gibbs) algorithm is as follows:   1. Specify the candidate function $q$ as required in the question. I multiplied each value in $(\boldsymbol{X}'\boldsymbol{X})^{-1}$ by $0.9$ for my specification of $\boldsymbol{\Omega}^{-1}$   2. Draw a candidate draw, $\boldsymbol{\beta}^c$, from $q(.)$   3. Calculate the probability $\alpha = \min\left(\frac{p(\boldsymbol{\beta}^c|\sigma^{(j-1)},\boldsymbol{y})/q(\boldsymbol{\beta}^c|\sigma^{(j-1)},\boldsymbol{y})}{p(\boldsymbol{\beta}^{(j-1)}|\sigma^{(j-1)},\boldsymbol{y})/q(\boldsymbol{\beta}^{(j-1)}|\sigma^{(j-1)},\boldsymbol{y})},1\right)$   4. With probability $\alpha$, set $\boldsymbol{\beta}^{(j)}=\boldsymbol{\beta}^c$, otherwise set $\boldsymbol{\beta}^{(j)}=\boldsymbol{\beta}^{(j-1)}$   5. Use $\boldsymbol{\beta}^{(j)}$ to define the conditional posterior of $\sigma$ from which $\sigma^{(j)}$ is drawn and so on. My code in R is:               ####################################################################################     ###########     Functions to be used in the MH algorithm      ######################     ####################################################################################          ## Specifying the conditional posterior pdf for beta (vector)          condbetavec = function(beta,sigma){            intcont = (2*pi)^(-3/2)*sigma^(-3)*(det(invxx))^(-1/2)              cond_pdf = intcont*exp(-(1/(2*sigma^2))*t(beta - betahat)%*%solve(invxx)%*%(beta - betahat))            return(cond_pdf)          }               ## Specifying the candidate pdf          cand = function(cand_beta,cand_sigma){            intcont = (2*pi)^(-3/2)*(cand_sigma)^(-3)*(det(cand_invxx))^(-1/2)              cand_pdf = intcont*exp(-(1/(2*(cand_sigma)^2))*t(cand_beta - betahat)%*%solve(cand_invxx)%*%(cand_beta - betahat))            return(cand_pdf)          }          ###########################################################     ########       Generation of data = y               #######     ###########################################################               true_beta = matrix(c(2,0.5,0.7),3,1)      ## True value of beta (vector) used in the data generating process (dgp)          true_sig = 0.3                            ## True value of sig used in the data generating process (dgp)               ## Setting the random number seed          set.seed(123456, kind = NULL, normal.kind = NULL)               ## Number of observations          nobs=20               ## Generating the values for x1, x2, x3 and y. Using Gaussian distribution here for convenience.          x1 = matrix(1,nobs,1)     x2 = matrix(c(rnorm(nobs,0,1)),nobs,1)     x3 = matrix(c(rnorm(nobs,0,1)),nobs,1)     u  = matrix(c(rnorm(nobs,0,1)),nobs,1)%*%true_sig          x = cbind(x1,x2,x3)          y = x%*%true_beta + u               ###########################################################     ### Specification of sample statistics     ###########################################################               invxx = solve(t(x)%*%x)     betahat = invxx%*%t(x)%*%y     cand_invxx = invxx*0.9  ## Specifying the var-covar matrix for the candidate pdf          df = nobs - 3     sighat = sqrt((t(y-x%*%betahat)%*%(y-x%*%betahat))/df)               ###############################################################     ###     ###  Generation of draws of beta/sigma/psi via Gibbs+MH sampling     ###     ###############################################################               B = 100+1        ## Extra +1 since the Gibbs sampler coded below starts at j=2, hence burn-in period is from j=2 to j=101     M = 1000         ## Number of draws after burn in period.     repl = B+M       ## Total number of Gibbs iterations including burn in and number of draws after burn in period.               ## Dimensioning the matrix of beta/sigma/psi draws ##                betav = matrix(0,repl,3)           sigv = matrix(0,repl,1)     psiv = matrix(0,repl,1)               ## Generating underlying normal random draws to be used in Gibbs+MH algorithm               stnbeta = matrix(c(rnorm((repl*3),0,1)),repl,3)                ## Specifying starting value for Gibbs+MH chain          sigv[1] = true_sig     betav[1,] = true_beta          ## Outer Gibbs Chain          for(j in 2:repl){            ## MH sub-step within the Gibbs to generate beta (vector) via candidate pdf            meanbeta = betahat       varbeta = (sigv[j-1])^2*cand_invxx            betagen = chol(varbeta)%*%stnbeta[j,] + meanbeta             if (runif(1,0,1)<min(((condbetavec(betagen,sigv[j-1]))/(cand(betagen,sigv[j-1])))/((condbetavec(betav[j-1,],sigv[j-1]))/(cand(betav[j-1,],sigv[j-1]))),1)){       betav[j,] = betagen       }else{       betav[j,] = betav[j-1,]       }                 ## Generation of sigma via its known inverted gamma conditional            a = t(y-x%*%betav[j,])%*%(y-x%*%betav[j,])            zvec = matrix(c(rnorm(nobs,0,1)),nobs,1)            chisq = t(zvec)%*%zvec            sigv[j] = sqrt(a/chisq)                 ## Draws of psi from its marginal posterior pdf            psiv[j] = (betav[j,2]+betav[j,3])/(sigv[j])^2                    }          #################################################################################     ## Estimation of the marginal posterior pdf of psi using kernel density smoothing     #################################################################################          ## Remove burn in period of psiv which is from j=1 to j=101 to create the vector psi which conists of 1000 draws after the burn in period           psi = psiv[102:repl]               ## Histogram of draws of psi          hist(psi,30,ylab="Percentage frequencies",xlab=expression(psi),          main=expression(paste("Histogram of draws of ",psi)))               ## kernel density plot of the draws of psi          d=density(psi)     plot(d,main=expression(paste("Kernel density estimate of marginal of ",psi)),          xlab=expression(psi),ylab=expression(paste("p(",psi,"|y)")))      Just wondering if anyone could verify whether the basic structure of my code is correct. Thanks.


title: Mixed-effects model with uneven samples

text: I have a somewhat trivial question regarding a mixed-effects model. The scenario is this:   1. 4 different treatments    2. 3 different time points   3. All time points have the same number of samples in each treatment group but the 3rd time point has half of the overall points ( _e.g._ , time points 1 and 2 have 40 points with 10 points per treatment group while time point 3 has 20 points with 5 points per treatment group). So what they are really interested in is whether there is a treatment effect. I thought I would model this as a nested ANOVA, fixed effect being treatment and time as a random effect (because presumably they could have taken readings at any time).               g.aov <- summary(aov(Y ~ treatment+Error(time/treatment), data=data))      I am concerned a bit about the sample imbalance in the third group and in addition I definitely see a time effect from the box plot so not sure if this model is capturing that well. Any advice would be appreciated.


title: How to check independence from empirical data

text: Let we've empirical data of two random variables X and Y. We can estimate empirically the degree of correlation between the two random variables by observing their scatter diagram, for instance. Since, **Independent random-variables $\Rightarrow$ zero-correlation** But **zero-correlation $\nRightarrow$ independence** **Question:** How to find the independence between two random-variables from empirical data?


title: Using a set of binary logistic regressions with multiple choice categorical response variable

text: I have categorical survey data on people's attitudes towards a certain policy area from 13 countries. The response variable is categorical, and includes 4 distinct answers that cannot be ordered. I would like to build a multi-level random-intercept and random-slope multinomial model. The problem is, that the number of level-2 cases is just 13, and the model does not converge, at least not in its multinomial form. So, as a second-best option, I am thinking about recoding the response variable into a binary form, run a series of multilevel logistic regressions, and then use predicted probabilities to show how the probability that a certain category of interest will be selected depends on my explanatory variables. This, apparently, is just a second-best option. I would like to know what the possible risks are of taking this approach, and what objections (from reviewers, supervisors etc.) should I expect.


title: How do I cluster documents using topic models?

text: Let us say I have a topic probability per document, for example:               Doc1: Science:0.6, History:0.2, Politics:0.1, Sports:0.1     Doc2: Science:0.3, History:0.5, Politics:0.1, Sports:0.1     Doc3: Science:0.8, History:0.1, Politics:0.05, Sports:0.05     Doc4: Science:0.2, History:0.2, Politics:0.4, Sports:0.2      It is fairly clear that Doc1 is similar to Doc3. If I use topic modeling as a dimensionality reduction tool, can I use clustering methods to cluster documents on topic space? What would be the appropriate clustering technique and distance function? Will K-means with Euclidean distance suffice? (It seems weird to calculate eucledian distance between two probability distribution) Also, topics might be correlated, what would be a way to handle that?


title: What is the correct way to test for significant differences between coefficients?

text: I'm hoping someone can help straighten out a point of confusion for me. Say I want to test whether 2 sets of regression coefficients are significantly different from each other, with the following set up:   * $y_i = \alpha + \beta x_i + \epsilon_i$, with 5 independent variables.   * 2 groups, with roughly equal sizes $n_1, n_2$ (though this may vary)   * Thousands of similar regression will be done simultaneously, so some kind of multiple hypothesis correction has to be done. One approach that was suggested to me is to use a Z-test: $Z = \frac{b_1 - b_2}{\sqrt(SEb_1^2 + SEb_2^2)}$ Another I've seen suggested on this board is to introduce a dummy variable for grouping and rewrite the model as: $y_i = \alpha + \beta x_i + \delta(x_ig_i) + \epsilon_i$, where $g$ is the grouping variable, coded as 0, 1. My question is, how are these two approaches different (e.g. different assumptions made, flexibility)? Is one more appropriate than the other? I suspect this is pretty basic, but any clarification would be greatly appreciated.


title: Coin flipping example

text: I have a bag of 10 coins; 1 of them has heads on both sides, while the rest are normal (have heads on one side and tails on the other). Each face of a coin is equally probable when tossed. You choose 1 of the coins from the bag at random and then flip it 100 times observing heads every time. What is the probability that you chose a normal coin? The probability that we take a normal coin is 9/10 and in 100 turns, it may be (9/10)^100, but I am not sure if this is correct.


title: GLM with Temporal Data

text: This is my first post on here, looking for some help. I am relatively new to analysis of temporal datasets. I have experience with R and developing linear models, so I am trying to figure out if the approach I am trying is appropriate with temporal data or if I need to do something different. Data: Monthly measurements of metal concentrations that span one year, with a few missing months. With the missing months, I am pretty I can't set up a true time series analysis in R (right?). I have several predictor variables (rainfall, temp, etc) that I want to evaluate. Approach: Construct a GLM with terms for each predictor and sampling month. I am considering using the gls function in nlme package in R so that I could examine autocorrelation. Goal: Identify factors to explain pattern in temporal water data Questions: - Can I use a glm approach with temporal data, even if there may be seasonality to the data? - How do i test for temporal autocorrelation in the error terms? I know there are some correlation structures that I could set up in gls but haven't figured out how to set this up and which one is best for testing temporal autocorrelation yet. Example dataset and approach is below. I would appreciate advice on if what I propose seems appropriate. If not, suggestions of alternative approaches would be a great help.               dat<-data.frame(month=c("Jan","Jan","Jan","Feb","Feb","Feb","Jun","Jun","Jun"),     wat=c(0.4,0.5,0.6,1,1.5,1.8,0.4,0.3,0.8),     rain = c(100,110,120,200,210,220,300,310,330),      temp = c(10,11,12,14,12,13,10,12,8))          m0<-gls(wat~month+rain+temp,correlation = corAR1(),data=dat)      Thanks!


title: How to sort a character vector numerically based on a trailing digit in R?

text: I have a vector `x = c("E1","E2","E4","E4","E5","E6","E10","I1","I2","I4","I12")` I need to sort that in a way the output should be `x = c("E1","I1","E2","I2","E4","E4","I4","E5","E6","E10","I12")` but when I use `sort(x)` it gives `c("E1","E10","E2","E4","E4","E5","E6","I1","I12","I2","I4")\` Is there anyway I could sort them as I wanted above?


title: Conditions for the M-estimator to converge to the true mean

text: Given i.i.d. samples from a gaussian distribution $X_1,...,X_n \sim N(\mu,\sigma) $ and the M-estimator, $\mu_m = \underset{a}{\operatorname{argmin}} \sum\rho(|X_i-a|)$, what properties on $\rho$ are sufficient to guarantee $\mu_m \rightarrow \mu$ in probability? Is $\rho$ being strictly convex and strictly increasing sufficient?


title: Stein's loss for multivariate normal covariance estimator

text: In James and Stein (1961) the authors consider the following loss function for an estimator $\hat{\Sigma}$ of the covariance matrix $\Sigma$ of a multivariate normal distribution: $$L(\hat{\Sigma}) = tr[\hat{\Sigma}\Sigma^{-1}] - \log|\hat{\Sigma}\Sigma^{-1}| - p.$$ This loss function has since been referred to as `Stein's loss' in several papers on regularized covariance estimation. Is there any intuitive justification for this loss function? I noticed that the loss function resembles the KL divergence between two multivariate normal distributions with the same means and with covariances $\hat{\Sigma}$ and $\Sigma$: $$2KL(N(0,\hat{\Sigma}) || N(0,\Sigma)) = \log|\hat{\Sigma}^{-1}\Sigma| - C\int_x [x^T (\hat{\Sigma}^{-1} - \Sigma^{-1})x ]\exp\\{-\frac{1}{2}x^T \hat{\Sigma}^{-1} x\\} dx$$ where $$C = (2\pi)^{-k/2}|\hat{\Sigma}^{-1}|^{1/2}$$ however I am not sure how to simplify the integral.


title: Siginificant size for a Control Group

text: I am designing project which requires me to select a Control Group. One the business requirements stipulates that I select as small a sample size as possible for the Control Group so there will be more records for the Test Group. I would like to select as small as possible a Control Grouo and still be confident that the result would be statistically significant. Does anyone know which test would be appropriate for this? Kim-Yen


title: Test to identify change in median in a time-series

text: the question is trivial, but is there any statistical test (can be done in R) to show that the median is changing in a time-series?? For example, if you go to the following link, you would notice that, after the first few samples (6 or 7), the time series takes a sharp rise. I was wondering if there is a formal statistical way of saying it! https://docs.google.com/file/d/0B6GUNg-8d30vMzB3cUxJM1dBbXc/edit?usp=sharing


title: Types of regression that might work on small samples

text: I know that this kind of discussion might be both a silly...very silly and a hot issue simultaneously...but I thought it would be nice to raise it : Which types of linear AND non-linear regression could be recommended when you have a small sample size (e.g. 50-70 records) and you want to predict a continuous outcome from 5-7 predictors (both categorical and continuous). Can k-fold cross-validated stepwise linear regression OR Partial Least Regression be a remedy to this problem? Are boosted trees OR Neural networks just insane even to think about them? Your responses are GREATLY WELCOMED and VERY MUCH APPRECIATED. THANK YOU!


title: Manipulating Binomial Distribution

text: Recently, I've been reading Yudi Pawitan's book, In All Likelihood. In the book, there's a section on profile likelihood; the methods explored in this section are subsequently applied to some data on heart attack prevalence amongst two distinct groups: people taking aspirin and people subjected to a placebo; the groups are modelled as $\text{Bin}(n_{a}, \theta_{a})$ and $\text{Bin}(n_{p}, \theta_{p})$ respectively. Since $n_{a}$ and $n_{p}$ are large while the event rates are small, $X_{a}$ and $X_{p}$ are approximately Poisson with parameters $n_{a}\theta_{a}$ and $n_{p}\theta_{p}$ respectively. You can see the example here (Pages 87 and 88). At the end, it is shown that the profile likelihood for a parameter of interest $\theta$ is given by: $$ L(\theta, \theta_{p}) = \text{constant} \times e^{-\theta_{p}(n_{a}\theta+n_{p})} \theta^{x_{a}} \theta_{p}^{x_{a}+x_{p}} (1)$$ where $\hat{\theta}_{p} = \frac{x_{a}+x_{p}}{n_{a}\theta+n_{p}}$, the MLE for $\theta_{p}$, should be substituted into the above model for $\theta_{p}$. and $n_{a}$ : number of people in aspirin group, $n_{p}$ : number of people placebo group, $x_{a}$ : number of heart attacks amongst aspirin group, $x_{p}$ : number of heart attacks amongst placebo group. also $\theta = \frac{\theta_{a}}{\theta_{p}}$, $\theta_{a}$ : probability of heart attack in aspirin group and $\theta_{a}$ : probability of heart attack in placebo group. The above likelihood is then further expressed as binomial, taking the following form: $$ L(\theta) = \text{constant} \times \left(\frac{n_{a}\theta}{n_{a}\theta+n_{p}}\right)^{x_{a}} \left(1-\frac{n_{a}\theta}{n_{a}\theta+n_{p}}\right)^{x_p} (2)$$ I have tried to understand how to get from the expression given in $(1)$ to the expression given in $(2)$, but when I make the substitution for $\hat{\theta}_{p}$, I get the following: $$ e^{-(x_{a}+x_{p})}\theta^{x_{a}}(x_{a}+x_{p})^{x_{a}+x_{p}}\left(\frac{1}{n_{a}\theta+np}\right)^{x_{a}+x_{p}} $$ where $e^{-(x_{a}+x_{p})} = \text{constant}$ as it doesn't depend on $\theta$. I believe this is along the right lines, but I'm not sure where to go next. Perhaps somebody could explain how to get from expression $(1)$ to expression $(2)$.


title: DNA use in court cases

text: I am currently studying the following case of Neil Owen, based on the following article I found a newspaper: > "A 20-year-old student was jailed for life yesterday for the brutal rape and > murder of a schoolgirl, after one of the biggest DNA testing programmes in > British criminal history. Neil Owen was arrested a year after the murder > when his genetic fingerprint was matched with DNA found at the scene, > following a mass DNA screening of 2000 men on the estate. He lived just 100 > yards from the victim's house. Laboratory tests revealed the chances of > anyone else being the killer were 1 in 160 million." Now first of all I am aware that there is an issue with prosecutors fallacy here. Because the 1 in 160 million is interpreted as P(innocence|matching blood type evidence) when it actually refers to P(matching blood type evidence|innocence). But my question refers to the reasoning of the defence. Counsel for the defence pointed out that there are about 30 million males in the United Kingdom, and argued that the correct probability that Owen was guilty is about **16/19** , not high enough to convict beyond a reasonable doubt. So my two questions are _1\. How do you think that the figure 16/19 was calculated? (I am sure the population of 30 million and the probability of 1 in 160 million was used?)_ _2\. What implicit assumptions were made, and how reasonable are they?_


title: PCA on correlation or covariance?

text: What are the main differences between performing Principal Components Analysis on a correlation and covariance matrix? Do they give the same results?


title: Why do we need to normalize data before analysis

text: I'm doing principal component analysis on my dataset and my professor told me that I should normalize the data before doing the analysis. My question is: "Why?" :) What would happen If I did PCA without normalization? Why do we normalize data in general? Could someone give clear and intuitive example which would demonstrate the consequences of not normalizing the data before analysis etc.? Thank you for any help =)


title: Categorical dependent variable in between subject design

text: My colleague has an experiment with a 2x2 all between subjects design with categorical variables predicting a three-level categorical dependent variable (choise between A, B, C/ coded as -1 0 1). Is this coding correct? Can a dependent variable have three levels? What kind of analysis can he best conduct?


title: Comparing two curves with different x-axis points - appropriate test?

text: I have two curves and I want to be able to calculate the probability of these curves coming from different distributions or another appropriate statistic. Each curve is fitted through the mean of clusters of data (x's and o's in diagram) at differing points on the x-axis, each data cluster is also non- parametrically distributed. Is there a statistical test that would be appropriate to tell me how likely it is that the data/curves come from different distributions i.e. one is significantly different than the other?                       +   -         x             |    o        ~             |    -        x~             |    o         ~             |     -        ~~             |     --        ~~             |      -         ~             |      --         ~~~             |       -            ~~~             |       -               ~~    x             |        -o               ~~~ x             |         -                  ~~              |         o-                  x  ~              |         o-                  x   ~ ~~             x             |           --                        ~~~ ~~~~~~~~~~~             |            ----        o                         x             |                - ----  o                         x             |                      ----- ----               o             |                        o       - ---- - -  - ---             |                                               o             |             +----------------------------------------------------------+      I've looked here Comparison of two curves but I believe my problem has some distinct differences


title: Is system identification and learning examples of parametric or nonparametric methods?

text: Difference between parametric and nonparametric statistical tests and What is the difference between distribution free statistics/methods and non-parametric statistics? is a good read for beginners. So, does system identification of linear models like Auto regressive, Moving AVerage, ARMA And learning of neural networks by estimating its weights fall under Parametric or non parametric methods?


title: Changing logistic regression's loss function

text: We're using logistic regression to predict events probability. Logistic regression tries to minimize the residual variance (sum of squared residuals). However, in our specific problem we would like to use a different loss function (for the errors in the logistic regression). Does anyone knows any implementation (or how to implement) such change? (preferably in R or Python) Thanks!


title: How to calculate the "unfairness" for repeated random selection?

text: Assume there are 16 people in a group, and each day 11 people are randomly selected from the group. Over a year, some people will be selected more times than other people. The question is, assume the person selected most often was selected $M$ times, and the person selected least often was selected $m$ times, what is the distribution (or the expected value) of the quantity $M/m$?


title: How can I compare the deviation from CSR in different point patterns?

text: I recently discovered the many tools of point patterns analysis and this is a very interesting field. I read a lot about how to look for deviations from complete spatial randomness (CSR) and I am able to analyze a point pattern at different scales using first- and second-order functions. However, I notice that I still don't know a technique to compare different point patterns to each other. For example, if I have two study areas A and B, I am able to tell that points within A and B are clustered, but **how can I tell that clustering within A is more intense than clustering within B?** I know about Moran'I as a global indicator of spatial autocorrelation, but is it suitable to study point patterns where no Z data is associated to the points? (maybe this is a very well knwon issue, but I couldn't find the keywords to help me find something revelant). Thank you for your help!


title: Interaction are not significant for model and coefficients, but main effect is significant

text: I am doing a social science research on behaviour. A purchase intention. I have intention to purchase as dependent variable, attitude as predictor and extraversion as moderator.               DV = intention to purchase     step 1     attitude (center)     extraversion (center)     step 2     Attitude X extraversion (center)      The model summary show significant only in model 1 with very little R-Square change (0.003), model 2 non-significant. The ANOVA show significant for both (0.000). The Coefficients show significant only for attitude. The extraversion only show non-significant. The interaction Attitude X extraversion show non- significant. But when I try to plot a High-Low graph. They are not parallel and show some cross together. What should I interpret this? Can I interpret the graph with Non-Significant interaction-term? Thank you very much.


title: Vectorized evaluation of a simple ANN on a mesh grid

text: I have written some code that implements a neural network which solves the `XOR` problem. I have a function `h` that receives 2-D points and the network matrices at each layer and computes the output of the net on that point. Now, what I want is to visually examine the decision boundary on the rectangle `(0,1) (0,1)`. To do that, my plan is to compute the output of the net for each of the points in that grid. I could just do two nested loops and calculate the output, but I'm looking (for some strange reason) for a more elegant vectorized way of doing that. I'm using `numpy` and `python`, but a `matlab` solution will be clear to me too. Any ideas? or am I fighting ghosts here?


title: Attribute values for minimizing outcome

text: I have data of the following format: Attr1,Attr2,Attr3....Attrn,Outcome (All attributes and outcome are numeric) I want to find out: 1\. With what values of each attribute will the Outcome be minimum? 2\. What attributes are most responsible for minimizing the Outcome?


title: Why is the expection of $E(Y_{T+1}|\Omega_T)$ greater than or equal to its previous value?

text: Consider the following model for $Y_t$: $\Delta$log($Y_{T+1})$ = $u_T$ where $u_T$ ~ IID Normal(0,$\sigma^2$). I want to forecast $Y_{T+1}$. Taking exponentials and then expectations, we see that $E(Y_{T+1}|\Omega_T)$ = $Y_T E(e^{u_{T+1}}|\Omega_T)$ = $Y_TE(e^{u_{T+1}})$ since $u_{T+1}$ is IID. Also, since $u_t$ ~ Normal(0,$\sigma^2$), we know (from the moment generating function) that $E(e^{u_t})$ = $e^{(1/2)\sigma^2}$ Question: Why is there a tendancy for the one step ahead forecast to be above its previous value, and not below. In mathematical terms, why is $E(Y_{T+1}|\Omega_T) \geq Y_T$?


title: Prediction using SVD and Fisher's linear discriminant

text: Where can I get an explanation of the procedure used when making a prediction using SVD? Let me elaborate a bit more. Suppose you have data in a matrix $A$ containing two classes. In particular, you could have $m$ attributes and $n$ objects. Half of the objects belong to the first class and the other half belong to the second class. Using SVD we obtain the matrices $U$, $S$, $V^{*}$. From the factor interpretation, we can regard the rows of matrix $U$ as a different view of the $m$ attributes of matrix $A$ and the columns of $V$ as a different view of the objects of matrix $A$. Therefore, it's possible to separate the objects depending on their class from $V^{*}$. For example, supposing that $V^{*}$ looks like this: $$ V^{*} = \begin{pmatrix} v_{1,1} & v_{1,2} & \cdots & v_{1,n} \\\\\ v_{2,1} & v_{2,2} & \cdots & v_{2,n} \\\\\ \vdots & \vdots & \ddots & \vdots \\\\\ v_{n,1} & v_{n,2} & \cdots & v_{n,n} \end{pmatrix} $$ Then we can separate the first $n/2$ columns and regard them as an alternative view of objects from the fist class. We can do the same with the last $n/2$ columns as an alternative view of objects from the second class. In practice, I have seen we take the matrix $\Sigma V^{*}$ of the decomposition: $$ \Sigma V^{*} = \begin{pmatrix} \sigma_{1} v_{1,1} & \sigma_{1}v_{1,2} & \cdots & \sigma_{n} v_{1,n} \\\\\ \sigma_{2} v_{2,1} & \sigma_{2}v_{2,2} & \cdots & \sigma_{2}v_{2,n} \\\\\ \vdots & \vdots & \ddots & \vdots \\\\\ \sigma_{n}v_{n,1} & \sigma_{n}v_{n,2} & \cdots & \sigma_{n}v_{n,n} \end{pmatrix} $$ Now that I think about it, I don't know why. Now we have enough information to calculate the mean of each class in this rather unusual view of our objects. So you can use Fisher's discriminant analysis to get projections of elements from each class in such a way to minimize the variances within elements of a given class and maximize the means between classes. This is succinctly described by the following equation: $$\vec{w} = \text{ arg max }_{\vec{w}} \displaystyle \frac{\vec{w}S_{b}\vec{w}}{\vec{w}S_{w}\vec{w}}$$ The solution is found by solving the associated generalized eigenvalue problem. This produces a set of eigenvectors and eigenvalues. Apparently, we need to choose the eigenvector corresponding to the largest eigenvalue. Why? After getting this $\vec{w}$ we have everything we need to start making predictions. For simplicity, consider a new object $\vec{b}$ that has the same kind of attributes as $A$. Supposedly, a prediction starts by multiplying $U^{T}$ and $\vec{b}$: $$ U^{T}\vec{b} = \begin{pmatrix} u_{1,1} & u_{1,2} & \cdots & u_{1,n} \\\\\ u_{2,1} & u_{2,2} & \cdots & u_{2,n} \\\\\ \vdots & \vdots & \ddots & \vdots \\\\\ u_{m,1} & u_{m,2} & \cdots & u_{m,n} \end{pmatrix} \begin{pmatrix} b_{1} \\\ b_{2} \\\ \vdots \\\ b_{n} \end{pmatrix} $$ Here the columns of $U^{T}$ represent attributes and the rows represent hidden factors. Likewise, the vector $\vec{b}$ represents a single object with attributes in the rows. As a result of this multiplication, we have an object described in terms of $m$ factors. Finally, to predict a class for our new object $\vec{b}$, we multiply $\vec{w}$ by the previous result: $$w^{T}U^{T}B = \begin{pmatrix} w_{1} & w_{2} & ... & w_{m} \end{pmatrix} \begin{pmatrix} u_{1,1} & u_{1,2} & \cdots & u_{1,n} \\\\\ u_{2,1} & u_{2,2} & \cdots & u_{2,n} \\\\\ \vdots & \vdots & \ddots & \vdots \\\\\ u_{m,1} & u_{m,2} & \cdots & u_{m,n} \end{pmatrix} \begin{pmatrix} b_{1} \\\ b_{2} \\\ \vdots \\\ b_{n} \end{pmatrix} $$ Since $\vec{w}$ is the vector that works as a projection to produce a decision threshold, I think it should live in the same space that other points from $V^{*}$. Therefore, its columns should represent factors. The result of $w^{T}U^{T}B$ produces a number. The comparison of this number with the threshold allows to make a decision on the predicted class. All of this seems rather reasonable. These operations work well and, of course, produce valid operations. However, I want to have a better grasp of this procedure.


title: Fitting a best fit line to this time series

text: I have the following hourly time series data and would like to fit a best fit line to it: ![enter image description here](http://i.stack.imgur.com/B03jZ.png) There seems to be a periodicity on a daily basis and a weekly basis. By this, I mean there are patterns that repeat every day (e.g. peaks during 7PM) and patterns that repeat every week (e.g.peaks during Sat, Sun). Here is a acf plot with a max lag of 168 observations (=hours): ![enter image description here](http://i.stack.imgur.com/8gLTX.png) There seems to be a lot of pattern here but I don't really know how to exploit it. I have tried using `auto.arima` from the `forecast` package:               > auto.arima(login.hourly)     Series: login.hourly      ARIMA(5,1,3)                              Coefficients:              ar1      ar2     ar3     ar4     ar5      ma1     ma2      ma3           2.3955  -2.0707  0.5538  0.0625  0.0028  -2.6221  2.5771  -0.9534     s.e.  0.0258   0.0521  0.0297  0.0160  0.0151   0.0119  0.0202   0.0092          sigma^2 estimated as 42.33:  log likelihood=-4815.75     AIC=9645.58   AICc=9645.7   BIC=9693.17     Warning message:     In auto.arima(login.hourly) :       Unable to fit final model using maximum likelihood. AIC value approximated      I would like to:   1. Fit and plot a best fit line along with the original data. I'm more interested in fitting a line to describe my data as opposed to forecast it.   2. Calculate some sort of statistic to tell me goodness of fit (RMSE or is there sth better?)


title: Confusion over Lagged Dependent and HAC Standard Errors

text: It seems that the consensus is that after fitting a linear regression via OLS, when you have included a lagged dependent variable (because it is supported by theory or appears needed based on the data) such as the following AR(1): $y_{t}=\beta_{0}+\beta_{1}y_{t-1}+\beta_{2}x_{1}+e_{t}$ the resulting coefficients are 1) Always biased 2) If the stationary condition holds ($|\beta_{1}|$<1) and there is no serial correlation remaining after the lagged DV is added then the coefficients are consistent. **Question 1: Is this summary correct?** **Question 2: If there is no serial correlation in the residuals after adding the lagged (1 or more lags) dependent variable do you typically need a standard error robust to serial correlation like Newey-West?** I would think the answer is no because there is no serial correlation but Wooldridge (page 431 of 3rd edition of introduction) points out that a serial correlation robust standard error "...is valid to use ....in models with lagged dependent variables..". Since I understand that Newey-West standard errors are only appropriate with consistent estimators, this must mean that after the lagged dependent is added, there is no more serial correlation. That makes me wonder why pointing out the ability to use the standard errors with models with lagged DV is relevant? **Question 3: As a general rule, is it better to just NOT use a lagged DV but instead use OLS to get non-biased estimates and then use a Newey-West or similar procedure for inference?** **ADD AFTER ALECOS RESPONSE:** Here is a snippet of the page ![enter image description here](http://i.stack.imgur.com/DCnuv.png)


title: Unbiasedness condition in ordinary kriging and simple kriging

text: I have this confusion. In ordinary kriging we have used the unbiasedness condition which gave the sum of weights equal to one. However, in the case of simple kriging we have no such conditions why? I used ordinary kriging and simple kriging in my dataset and I have a huge difference in kriging in predicting the values. AFAIK, simple kriging is a linear regression which ordinary kriging is a constrained linear regression where sum of weights equal to 1.


title: Is there a PDF for a generalized non-central chi-squared distribution

text: Is there a PDF for a distribution defined as a sum of squares of random variables pulled from a family of normal distributions with different standard deviation? Is there a way of analytically writing the PDF of such a distribution as a hypergeometric function without prior normalization to make all the underlying normal distributions have unit variance? A possible case of application would be an N-dimensional length of displacement from center under the effect of a random, normally distributed forces along each dimension, each with it's own characteristic variance and mean, different between dimensions.


title: Distribution of sum of squared normal distributions with different means and variances?

text: Let $ Z=X_1^2+X_2^2\cdots\cdots X_j^2 $, such that $X_i \sim \mathcal{N}(\mu_i,\sigma_i^2)$. All $X_i$'s are independent of each other and $\mu_i$ and $\sigma_i^2$ are all different from one another. Question is: 1) What is the distribution of Z ? 2) If it is non-central Chi square distribution, what it ll be in terms of means and variances? Note: I'm not going to call $X^2$ or $Y^2$ as Chi Square distribution, as that complicates the matter. Since in Central $\mu_i=0$ and in non-central Chi square distribution each term in divided by respective variance, that don't represent our case.


title: Finding the experimental region of a linear model

text: Given a set of data with many independent variables and a single dependent variable how can I find the "experimental region" of the resulting model? This is the region within which we can make valid predictions and doesn't go outside of the bounds of the data we were provided. Also what is a good way to articulate this region?


title: compare binary responses with ties

text: Suppose I am testing two drugs, I and II. The response profile is either I > II, I < II or I = II. I was trying to figure out if I is in fact better than II. One way I thought of was to throw out all the ties (I == II), and then assign a 1 where I > II and a 0 where I < II. Then this boils down to a binomial distribution (which I can approximate with the normal). But I was wondering if throwing away the ties is really appropriate... If I don't throw away ties, I can't use the binomial distribution (so would I switch to some multinomial test)? The other thing I can think of is some wilcoxon rank sum test or mann whitney u test with ties... but I was wondering if anyone else has an idea? The concern here is that there will be a LOT of ties... Thanks!


title: permutation/bootstrapping for testing difference of difference

text: I have a situation where I have 3 treatments/conditions A, B, C. For each condition I have multiple trials. I compare the means of trials under each condition A, B and C. This gives me 3 values. Then I am computing a difference A-B and A-C. I need to obtain confidence interval/error bars for the single value obtained when I do the subtraction of mean response to condition A and mean response to condition B (A-B) same goes for A-C. Essentially I want to ask is the difference seen in mean response A-B different from mean response A-C. I know potentially some bootstrapping/permutation test is what is needed but not sure which way to do this since its not a direct comparison of mean of A vs mean of B but rather A-B compared to A-C. Help will be appreciated.


title: How to use Viterbi algorithm for multiple chain HMM?

text: I need to build a HMM with two chains: the structure is as the figure below. How can I use the Viterbi algorithm for this model? ![enter image description here](http://i.stack.imgur.com/mW1gB.jpg)


title: Maximum entropy and non-informative distribution

text: Is Maximum Entropy rule equivalent to non-informativeness? In other words, when maximizing the entropy of a distribution, given some known stuff, is it equivalent to finding to most non-informative distributions? If so, how do you justify this mathematically?


title: How to test which data match model at best

text: ### Description I have 1D data with $N$ normally distributed clusters. I have to find a cluster, which is the worst (differs at most from the normal distribution). ### My approach I calculate $sq = \frac{(f(x) - y)^2}{\\#y}$, where $f(x)$ is value of normal PDF with mean equal to the center of the cluster and sigma equal to cluster's "radius", $\\#y$ is the total number of observations. Cluster with highest value of $sq$ is considered to be the worst one. ### Question Problem is that number of points per cluster differs a lot (one cluster could have 3000 points and other 300). And imho I think that if I had small errors and many points, I would end up with larger $sq$, then if I had bigger errors and small amount of points. Can you point me the right way?


title: How do I estimate the mean of a variable for different groups by controlling for age, gender, education...?

text: I am trying to estimate the mean of a variable for 2 different groups. If in Stata I use               . bys group: sum variable       I'll get the mean. However, there are differences among two groups in terms of age, gender, education... And I have to control for that. To be more clear, let's say my groups are immigrants and natives. I can get the mean by using               . sum variable if immigrants==1      . sum variable if immigrant==0       However, the characteristics of the two groups are different, i.e. immigrants are younger, less educated, etc. Hence, the difference between simple means across groups can be a result of being a immigrant but can also be because of different characteristics. Therefore, I have to control for the group characteristics among natives and immigrants. One way to do that is using matching immigrants with natives with similar characteristics. Then, I calculate the mean for immigrants and similar natives. I was wondering if there is another way of doing this without using any matching methods but with some simple Stata commands.


title: How do difference-in-difference designs account for temporal autocorrelation

text: Although there are doubtless many techniques for studying the impact of a discrete intervention over time, I am interested in two which have achieved widespread adoption in the social sciences:   * Interrupted time series   * Difference-in-difference The former applies all the principles of time series analysis, using ARIMA models to account for non-stationarity, autocorrelation, etc. The latter uses linear regression or variants, but allows for a control group. If anything, difference-in-difference is more widely adopted, particularly in the economic literature (after Card, Angrist, etc.) but also in health-related fields and possibly education (where it seems to be called "comparative interrupted time series"). Yet in the presence of substantial autocorrelation (likely, in the types of datasets being used), it would seem to over-state the precision with which results are estimated. My questions are thus:   1. How are difference-in-difference analyses valid if autocorrelation exists?   2. Is there a time-series technique that allows a control group to be utilized?


title: Confused with basic time series terminology

text: I am a bit confused with the basic time series terminology: Consider the following words:   * fitted values   * forecasted values   * in-sample forecasts   * out-of-sample forecasts   * in-sample fit I am using GARCH models in my thesis. So the equations have a lagged character. E.g. GARCH(1,1): $\sigma_t^2=\alpha_0+\alpha_1\epsilon^2_{t-1}+\beta_1\sigma_{t-1}^2$ My fitted values are the $\hat{\sigma}^2_{T+1|T}$ with using in-sample data, i.e. I have estimates of the alpha and beta (estimation using all in-sample data) and I plug in the in-sample data. Due to the lagged character of the equation this could also be called in-sample forecasts, right? This would be $\hat{\sigma}^2_{T+1|T}$, $\hat{\sigma}^2_{T+2|T}$, $\hat{\sigma}^2_{T+3|T}$, $\hat{\sigma}^2_{T+4|T}$.... Forecasted values in general means forecasting the volatility equation at the end of the time horizon, so this is out-of-sample, right? This would be out- of-sample forecasts? I am not sure about the meaning "in-sample fit". I mean, fit is always in- sample right? Or otherwise this would imply that there is an out-of-sample fit, which makes no sense? Does the word combination in-sample fit exists?? Is it right, that I can use fitted and forecasted values or in-sample forecasts and out-of sample forecasts? Because I think that fitted values is the same as in-sample forecasts and forecasted values is equal to out-of- sample forecasts


title: How to design a 2-level fractional factorial experiment given some prior knowledge

text: We're designing a trial to test four factors (call them `A`, `B`, `C` and `D`) for an experiment. Each factor has two levels: `+` and `-`. The `+` of factor `A` is believed to lead to better results and the `-` of factor `B` is also believed to lead to better results. We know nothing about `C` and `D`, nor whether any interaction effects exist or not. We want to get the following goals:   1. Find the combination of factors to get the best result.   2. Build a model containing `A`, `B`, `C` and `D` to predict the result of experiments. (This requires that `A` and `B` should still be in this design) How can we get a fractional factorial design based on our prior knowledge of the effects of `A` and `B`?


title: Measuring NN saturation: calculating probabilities

text: I am trying to devise a measurement of neural network saturation for my NN saturation study. Some background: saturation occurs when a hidden neuron outputs values close to the extremes (usually 0 and 1), and practically nothing in-between. Neuron paralysis occurs when a neuron always outputs a similar value, regardless of the input. Neurons themselves are simple sigmoidal functions that take input (any real number, a result of a weighted sum of many inputs) and output a single value in a bounded range (typically [0,1]). I would like to be able to somehow measure the extent of saturation and paralysis. Can statistics be applied? Given the hidden neuron outputs for all patterns in the training set, would it be possible (and legit) to calculate the probability of a hidden unit outputting a value close to the boundary, or the probability of a neuron outputting the same value all the time? I suppose every pattern (input) can be treated as an independent event, although the neuron output will obviously depend on the weights applied to the inputs the sum of which is fed into the neuron - would this make it a case of conditional probability?


title: Calculating means and confidence intervals for groups with multiple observations in each subject using the NLME package in R

text: I need your help to identify group means and calculate the corresponding 95% confidence intervals for a set of independent samples with 2 dependent observations within each of those samples. I'm using R with the `NLME` package. Below is a simplified example with a similar structure to my actual data. Let's assume that both wheels of a total of 7 independent bicycles (n=7) were either equipped with very robust tires (n=4 with `tire=1`) or standard tires (n=3 with `tire=0`), both with a profile of 10.0mm depth. All bikes were used to ride the same cumulative distance, but based on the riding style and road conditions, the tires mounted on one bike might be more similar to each other, than the ones on a different bike. After the distance was reached, the profile depth was measured (possible range 0.0 to 10.0 mm) on all tires (a total of 14 observations). Below is the described dataset:               > example        ID tire depth     1   A    0  3.67     2   A    0  3.67     3   B    0  1.33     4   B    0  3.33     5   C    0  1.33     6   C    0  1.67     7   D    1  6.17     8   D    1  4.00     9   E    1  2.37     10  E    1  5.50     11  F    1  3.50     12  F    1  4.67     13  G    1  5.17     14  G    1  5.83      In order to assess whether `tire 1` is better than `tire 0` in terms of wear, I've used the `lme` function of the `NLME` package.               > summary(lme(depth ~ tire, data = example, random = ~1 | ID))      I've interpreted the following output as a mean difference of 2.15mm in favour of the robust tire (p=0.027).               Linear mixed-effects model fit by REML      Data: example             AIC      BIC    logLik       51.09194 53.03156 -21.54597          Random effects:      Formula: ~1 | ID             (Intercept) Residual     StdDev:   0.3138429 1.207888          Fixed effects: depth ~ tire                    Value Std.Error DF  t-value p-value     (Intercept) 2.50000 0.5253550  7 4.758687  0.0021     tire        2.15125 0.6949794  5 3.095415  0.0270      Correlation:           (Intr)     tire -0.756          Standardized Within-Group Residuals:            Min         Q1        Med         Q3        Max      -1.8180875 -0.7997771  0.2085848  0.8333573  1.2146423           Number of Observations: 14     Number of Groups: 7       Using the DF and the Std. Error of the output I'm able to calculate the 95% CI of the effect of the robust tire when compared with the standard tire (i.e. 2.15mm, 0.36 to 3.94mm). However, I would like to describe the means and the corresponding 95% CI of each individual group. That is, we observed a profile depth of 4.65mm (X.XX to X.XX) in the robust tires and 2.5mm (X.XX to X.XX) in the standard tires, in a way that still takes into account that not all observations are independent.


title: What does "curvilinear" mean?

text: As far as I can tell, _curvilinear_ is defined vaguely but means the same as _nonlinear_. Is that correct? Or does _curvilinear_ have a distinct definition?


title: Derivative of Marginal df equal to marginal pdf?

text: I am curious as to when the derivative of a marginal df is equal to the marginal pdf for a (continuous) random variable. As a simple example, let Z = (X,Y) be a continuous random vector and F be the df of Z. Assuming the processes of integrating and differentiation can be exchanged, then it would seem like the marginal pdf for X should be equal to the derivative of the marginal df for X. Is this correct and if so, when can we say this is possible (I believe that for 1 variable, F uniformly continuous generally implies derivatives and integrals can be exchanged with limiting processes). Thank you, Matt


title: mixed between-within design when covariate-DV relationship is non-linear (quadratic)

text: My question is about what method to use for a mixed between-within design when the relationship between the covariate and the DV is non-linear (quadratic). My DV is reasoning task score. There are two types of questions in the task - let's call them "control" and "test" questions. My IVs are question type, age (which ranges from 11-50 years, but with uneven frequency distribution - more adolescents) and gender. N is roughly 170. I'm interested in how age - especially among the adolescents - moderates performance on the reasoning task and how this might differ between the two question types. I anticipate that performance on the control questions would increase more rapidly with age than for the test questions (i.e., interaction between age and question type). I'm also expecting that females will be better than males on the test questions (and possibly the control questions too) but am not sure whether or how gender might interact with age. A mixed design ANOVA without age (i.e., just question type and gender) shows significant main effects of both factors but no interaction. I had planned to use a mixed design ANCOVA with gender as a between-subjects factor, age as a covariate (mean centred) and question type as a repeated measure. However, the intial difficulty that I'm having is that the scatterplots suggest an inverted U shape relationship between age and performance on the task (for both question types). Obviously this means it's inappropriate to proceed with a model that assumes linear relationships between the covariate and the DV. So here's my question: is there a way to approach this using a nonlinear/nonparametric method? I haven't been able to find anything that suggests a nonparametric method for a mixed between-within design. I'd prefer something that can be implemented in SPSS if possible as that's the package I'm familiar with - but I am willing to consider other options if they are accessible to a newbie.